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1  Introduction 


Screen/film  mammography  is  widely  recognized  as  being  the  only  effective  imaging 
modality  for  the  early  detection  of  breast  cancer  in  asymptomatic  women  [1].  Screening 
asymptomatic  women  using  screen/film  mammography  has  been  shown  to  signihcantly 
reduce  breast  cancer  mortality  [2].  Breast  cancer  currently  accounts  for  32%  of  cancer 
incidence  and  18%  of  cancer  mortality  for  women  in  the  United  States.  There  were  182,000 
new  cases  of  breast  cancer  in  the  United  States  in  1993  and  46,000  deaths.  Five  year 
survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being  localized, 
falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [3].  The  early  detection 
of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce  breast  cancer 
mortality. 

Major  advances  in  screen/him  mammography  have  occurred  over  the  past  decade  [4] 
which  have  resulted  in  significant  improvements  in  image  resolution  and  film  contrast.  Of 
major  importance  is  that  these  improvements  have  been  achieved  at  reduced  radiation 
doses.  Despite  these  advances,  however,  screen/film  mammography  remains  a  diagnostic 
imaging  modality  where  image  interpretation  remains  very  difficult.  Breast  radiographs  are 
generally  examined  for  the  presence  of  malignant  masses  and  indirect  signs  of  malignancy 
such  as  the  presence  of  microcalcifications  and  skin  thickening.  Unfortunately,  it  is  unlikely 
that  major  improvements  in  imaging  performance  will  be  achieved  by  technical  advances  in 
screen/film  radiography  alone. 

The  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  minor 
difference  in  x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [5]. 
This  fact  makes  the  detection  of  small  malignancies  problematical,  especially  in  younger 
women  who  have  denser  breast  tissue.  Although  calcifications  have  high  inherent 
attenuation  properties,  their  small  size  also  results  in  a  low  subject  contrast  [6].  As  a 
result,  the  visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be 
a  problem  in  mammography  as  it  is  currently  performed  using  analog  film. 

Improvements  in  the  ability  of  screen/film  mammography  to  detect  small  tumors  and 
microcalcifications  is  more  likely  to  occur  by  improving  the  visibility  of  these  features.  It 
has  been  suggested  that  as  normally  viewed,  mammograms  display  only  about  3%  of  the 
information  they  detect!  [7]. 

Our  approach  to  feature  analysis  and  classification  is  motivated  in  part  by  recently 
discovered  biological  mechanisms  of  the  human  visual  system  [8].  Both  multiorientation 
and  multiresolution  are  known  features  of  the  humair  visual  system.  There  exist  cortical 
neurons  which  respond  specifically  to  stimuli  within  certain  orientations  and  frequencies. 
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In  this  report  we  describe  exciting  new  results  accomplished  during  our  second  year  of 
study.  In  addition,  we  have  continued  our  efforts  in  the  development  of  wavelet  transforms 
that  exploit  orientation  and  frequency  selectivity  to  make  mammographic  features  more 
obvious  through  localized  contrast  gain.  Below,  we  describe  in  executive  summary,  recent 
accomplishments  made  in  four  subtask  areas  of  our  research  and  development  project.  The 
methodology,  approach  and  experimental  studies  are  reported  in  detail  in  the  body  of  this 
report. 

1.1  Overview  of  Contents 

Contrast  enhancement  by  multiscale  nonlinear  operators. 

We  first  present  a  systematic  approach  to  accomplish  image  contrast  enhancement  by 
multiresolution  representations  of  the  dyadic  wa,velet  transform.  This  work  expands  upon 
our  initial  algorithms  development  during  the  first  year  of  our  investigation.  We  formulate 
two  examples  in  which  a  linear  enhancement  operator  is  shown  mathematically  equivalent 
to  traditional  unsharp  masking  with  a  Gaussian  low-pass  filter.  A  formal  analysis  of  wavelet 
filter  selection  and  associated  artifacts  is  carried  out.  We  show  that  transform  coefficients, 
modified  within  each  level  of  scale  by  nonlinear  operators,  can  make  more  obvious  unseen 
or  barely  seen  features.  In  addition,  we  introduce  an  edge-preserved  denoising  stage  based 
on  wavelet  shrinkage  with  adaptive  thresholding,  and  demonstrate  that  noise  suppression 
and  contrast  enhancement  can  be  achieved  simultaneously  within  the  same  framework. 

An  interactive  tool  for  editing  digital  mammograms. 

The  second  section  of  our  report  describes  the  design  and  functionality  of  XMam,  a 
computer  image  editor  that  allows  radiologists  to  indicate  on  digitized  mammograms 
regions  considered  as  having  signs  of  cancer.  This  is  accomplished  interactively  on  the 
computer  screen  by  enclosing  suspicious  regions  on  the  digitized  images  with  polygons. 

The  purpose  of  this  tool  is  to  develop  a  database  of  cases  diagnosed  with  cancer  based 
on  mammographic  screening,  and  later  confirmed  or  denied  as  positive  by  biopsy.  This 
shall  allow  us  to  compile  a  large  number  of  positive,  borderline  (difficult  to  diagnose)  and 
negative  cases  and  provide  ground  truth  for  the  development  and  verification  of  our 
computer  algorithms  for  the  detection  of  breast  cancer. 

Local  feature  analysis  via  interval  wavelets. 

This  section  of  the  report  introduces  a  novel  approach  for  accomplishing  interactive 
feature  analysis  by  overcomplete  multiresolution  representations.  Traditional  wavelets 
adapted  to  “life  on  an  interval”,  can  overcome  “edge  effects”  of  wavelet  representations  on 
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a  line.  Methods  of  contrast  enhancement  are  described  based  on  two  overcomplete 
multiresolution  representations  (interval  wavelets):  (1)  Deslauriers-Dubuc  interpolation,  (2) 
Average  interpolation. 

We  show  cjuantitatively  that  transform  coefficients,  modified  by  an  adaptive  non-linear 
operator,  can  make  more  obvious  unseen  or  barely  seen  features  of  mammography  without 
requiring  additional  radiation. 

Arbitrary  regions  of  interest  (ROl)  of  a  mammogram  are  enhanced  by  average 
interpolation  (Al)  and  Deslauriers-Dubuc  interpolation  (DD)  representations  on  an 
interval.  The  results  of  local  (ROl)  enhancement  and  global  enhancement  of 
mamniographic  features  are  compared  quantitatively.  We  demonstrate  that  our  method 
can  provide  radiologists  with  an  interactive  capability  to  support  high-speed  localized 
processing  of  selected  (suspicious)  areas  (lesions). 

These  results  augment  our  preliminary  findings  accomplished  during  the  first  year  of 
our  study,  and  further  demonstrate  that  features  extracted  from  multiscale  representations 
can  provide  an  adaptive  mechanism  for  accomplishing  local  contrast  enhancement.  In 
addition,  we  show  that  interval  wavelet  processing  can  be  carried  out  in  real-time  (at  video 
frame  rates)  over  selected  areas  of  arbitrary  shape.  This  is  significant  in  consideration  of 
the  large  image  matrix  sizes  produced  by  digital  detectors.  By  improving  the  visualization 
of  breast  pathology  we  can  improve  chances  of  early  detection  while  requiring  less  time  to 
evaluate  mammograms  for  most  patients. 

Quantitative  evaluation  of  clinical  images. 

The  final  section  of  our  report  relates  to  the  development  of  objective  ways  to  assess 
the  performance  of  wavelet  image  processing  algorithms.  The  objective  is  to  develop 
techniques  to  evaluate  wavelet  algorithms  so  they  can  then  be  optimized  for  clinical  use  in 
mammography. 

Substantial  progress  has  been  made  in  the  development  of  techniques  to  assist  in  the 
quantitative  evaluation  of  wavelet  based  image  processing  algorithms.  In  addition,  these 
techniques  have  been  applied  to  optimize  the  three  parameters  available  in  current  wavelet 
based  algorithms.  Specific  achievements  in  the  past  year  include: 

•  a  demonstration  of  the  ability  of  optimized  wavelet  based  algorithms  to  make  visible 
simple  objects  in  a  noisy  background  which  were  previously  invisible; 

•  a  demonstration  of  the  inherent  superiority  of  wavelet  based  algorithms  for  the 
detection  of  simple  objects  as  compared  to  algorithms  frequently  used  in  medical 
imaging  including  unsharp  mask  enhancement  and  median  filtering; 


6 


•  a  demonstration  of  the  special  requirements  of  wavelet  algorithms  when  enhancing 
the  visibility  of  features  of  specific  interest  in  mammography,  namely 
microcalcifications,  masses  and  fibril  structures. 

In  the  next  section,  we  shall  describe  in  detail  our  wavelet  processing  algorithms, 
experimental  methods  and  example  results  obtained.  In  addition,  we  list  and  summarize 
publications  and  presentations  made  by  our  researchers  during  the  past  year  of  our 
investigation. 
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Body 

2.1  Contrast  Enhancement  by  Multiscale  and  Nonlinear 
Operators 

2.1.1  Introduction 

Image  enhancement  techniques  have  been  widely  used  in  radiology,  where  subjective 
quality  of  images  is  important  for  reliable  diagnosis.  Image  contrast  is  a  important  factor 
in  subjective  quality.  Many  algorithms  for  contrast  enhancement  have  been  proposed.  A 
comprehensive  survey  is  presented  in  reference  [9].  Among  them,  histogram  modification 
and  edge  enhancement  techniques  are  most  commonly  used. 

Histogram  modification  techniques  [10,  11]  are  attractive  due  to  their  simplicity  and 
speed,  and  have  achieved  acceptable  results  for  some  applications.  In  general,  a 
transformation  function  is  derived  from  a  desired  histogram  and  the  histogram  of  an  input 
image.  Note  that  the  transformation  function  is  in  general  nonlinear.  For  continuous 
functions,  an  information  lossless  transformation  may  be  achieved.  However,  for  digital 
images  with  a  finite  number  of  gray  levels,  such  a  transformation  results  in  information  loss, 
due  to  quantization  errors.  For  example,  a  subtle  edge  may  be  merged  with  its  neighboring 
pixels  and  disappear.  Attempts  to  incorporate  local  context  into  the  transformation 
process  have  not  been  fruitful.  For  example,  simple  adaptive  histogram  equalization  [12] 
with  a  fixed  contextual  region  (window)  cannot  adapt  to  features  of  distinct  sizes. 

Most  edge  enhancement  algorithms  share  a  common  strategy  implicitly:  edge  detection 
and  subsequent  “crispening”.  “Unsharp  masking”  sharpens  edges  by  subtracting  a  portion 
of  a  Laplacian  filtered  component  from  an  original  image.  This  technique  was  justified  as 
an  approximation  of  a  deblurring  process  in  [13].  Loo  et  al.  [14]  studied  a  extension  of  this 
technique  in  the  context  of  radiographs.  Another  extension  based  on  Laplacian  filtering 
was  proposed  in  [15].  However,  these  (unsharp  masking)  techniques  remain  limited  by  their 
linear  and  single  scale  properties,  and  are  less  effective  for  images  containing  wide  range  of 
features  such  as  mammography.  In  an  attempt  to  overcome  these  limitations,  a  local 
contrast  measure  and  nonlinear  transform  functions  were  introduced  in  [16],  and 
subsequently  refined  in  [17].  However,  limitations  remain  in  these  nonlinear  methods:  (1) 
They  operate  on  a  single  scale,  (2)  There  is  no  explicit  noise  suppression  stage  (actually 
noise  may  be  amplified  significantly),  and  (3)  Nonlinear  transform  functions  were 
introduced  ad-hoc  without  a  rigorous  mathematical  analysis  of  their  enhancement 
mechanisms. 
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Figure  1:  One  dimensional  discrete  dyadic  wavelet  transform  (three-level  shown). 

Recent  advancement  of  wavelet  theory  has  sparked  researchers  in  the  applications  of 
image  contrast  enhancement  [18,  19,  20,  21,  22,  23,  24].  These  early  studies  showed 
promise,  but  remained  at  an  experimental  level.  In  this  report,  we  give  a  detailed 
mathematical  analysis  of  a  dyadic  wavelet  transform,  and  reveal  its  connection  to 
traditional  technique  of  unsharp  masking.  In  addition,  we  propose  a  simple  nonlinear 
enhancement  function  and  analyze  the  problem  of  artifacts.  Moreover,  we  introduce  an 
explicit  denoising  stage  that  preserves  edges  using  wavelet  shrinkage  [25]  with  adaptive 
thresholding. 

These  techniques  are  described  in  the  following  sections:  Section  2.1.2  presents  a  one 
dimensional  dyadic  wavelet  transform.  Section  2.1.3  analyzes  linear  enhancement  and  its 
mathematical  connection  to  traditional  unsharp  masking.  Section  2.1.4  analyzes  a  simple 
nonlinear  enhancement  by  point-wise  functional  mapping.  Section  2.1.5  introduces 
denoising  with  wavelet  shrinkage  and  an  adaptive  approach  for  finding  threshold  values. 
Section  2.1.6  presents  a  two  dimensional  extension  for  digital  mammography  and  special 
procedures  developed  for  denoising  and  enhancement  that  avoid  orientation  distortions. 
Section  2.1.7  presents  sample  experimental  results  and  comparisons.  Finally,  Section  2.1.8 
concludes  our  discussion  and  proposes  future  research  directions. 

2.1.2  One  dimensional  discrete  dyadic  wavelet  transform 

A  fast  algorithm  [26]  for  computing  1-D  discrete  dyadic  wavelet  transform  (DDWT)  is 
shown  in  Figure  1.  The  left  side  shows  the  structure  of  decomposition,  and  the  right, 
reconstruction.  For  an  N-channel  structure,  there  are  N  —  1  high-pass  or  band-pass 
channels  and  a  low-pass  channel.  Thus,  the  decomposition  of  a  signal  produces  —  1  sets 
of  wavelet  coefhcients  and  a  set  of  coarse  signal. 

Since  a  dyadic  wavelet  is  not  an  orthogonal  basis,  there  is  no  down-sampling  and 
up-sampling  shown  in  Figure  1,  and  its  decomposition  hlters  differ  from  its  reconstruction 
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Figure  2:  An  equivalent  multi-channel  structure  for  three-level  DDWT. 


filters. 

For  simplicity  of  analysis,  an  equivalent  multi-channel  structure  is  shown  in  Figure  2. 
This  computational  structure  also  makes  obvious  the  potential  for  parallel  processing. 

We  shall  refer  to  filters  and  in  Figure  2  as  forward  filters  and  inverse 

filters,  respectively.  Their  relationship  to  hlters  G(in),  K{u)  and  F[[u)  are  explicitly  given 

by 


N-\ 


Fo(oj)  =  Giio),  Fm{u;)  =  n 


1=0 


Fmi<^)  = 


'm—1 


n 


L  ;=o 


G{2^u;),  l<m<N-l. 


and 


N-l 

/o(cu)  =  Kiu),  iNiio)  =  n 

1=0 

rm—1 


I'm  (^) 


n  hf*(2'a;) 

L  1=0 


/i(2™a;),  1  <  m  <  iV  -  1. 


Since  filters  H{u)^  G(co)  and  K(u;)  satisfy  condition 


\\H{u;)f  +  Giio)K{u)  =  G 


filters  Fm{>^)  and  /^(cu)  completely  cover  the  frequency  domain, 


Channel  frequency  responses  Cm{ij)  can  be  written  as 


f  1  -  WHHf 

,  ?Ti  =  0, 

Cmi^)  — 

\  T-rm-l 

lb=o 

H{2^u) 

'  [l  -  ||i/(2”<.,)||t 

,  1  <  m  <  (iV  —  1), 

-i-rAf-l 

li(=0 

H{2^u) 

2 

^  m  =  N. 
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Figure  3:  Channel  frequency  responses  for  fV  =  6,n  =  1  and  (a)  p  =  0  and  (b)  p  =  1. 


As  an  example,  we  consider  an  extension  of  the  class  of  filters  proposed  in  [26] 


H{u)  = 


cos  — 


where  p  =  0,  or  1.  Let 


then  we  can  show  that 


m— 1 


JJ  cos(2^ 
l;=o 


sin(2”^-^u;)' 

L  2™sin(f)  '  ’ 


and  therefore 

0m,4n+2p(^)  Cm+l,4”+2p(^)  >  0  ^  fU  ^  (A^  1), 

0W,4n+2p(^)  •>  ^ 

Note  that  0o,n(<^)  =  1,  and  for  0  <  m  <  A", 


(1) 


(2) 


(3) 


,  X  2n+p— 1 

Cm{^)  =  0m,4n+2p(^)  ~  0!7i+l,4n+2p(a^)  =  sin  fw  j  4  0m,4n+2p+2  (<^)  [cOS  ^2  (jjj 

^  1=0 


21 


and  sin"  (l)  is  the  frequency  response  of  the  discrete  Laplacian  operator  of  impulse 
response  {1,  —2, 1}. 

Qm,qi^)  with  even  exponential  q  is  approximate  a  Gaussian  function,  while  the 
frequency  responses  of  channels  (0  <  ?u.  <  N)  are  approximate  a  Laplacian  of  Gaussian. 
Figure  3  shows  the  distinct  channel  frequency  responses,  and  Figure  4  compares  02,4(<^) 
and  02, 6(^)  with  related  Gaussians. 
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Figure  4:  (a)  02,4 (<^)  compared  with  the  Gaussian  function  e 
(h)  02,6(^)  compared  with  the  Gaussian  fnnction  . 


2.1.3  Linear  enhancement  and  unsharp  masking. 

Review  of  unsharp  masking 

A  protot5^pe  of  unsharp  masking  was  defined  [13]  as 

5(®,  y)  =  s{x,  y)  -  kAs{x,  y),  (4) 

where  A  =  ^  ^  is  the  Laplacian  operator.  However,  this  original  formula  worked  only 

at  the  level  of  finest  resolution.  More  versatile  formulas  were  later  developed  in  two  ways, 
described  below. 

One  way  to  extend  the  original  formula  is  based  on  exploiting  the  averaging  concept 
behind  the  Laplacian  operator.  The  discrete  form  of  the  Laplacian  operator  may  be 
written  as 

As(z,  j)  =  [s(f  +  l,i)  -  2s{iJ)  +  s{i  -  l,j)]  +  [^(f,  j  +  1)  -  2s{iJ)  +  s{i,j  -  1)] 

=  -5  |5(z,  j)  -  i  [s(f  +  l,i)  +  5(f  -  1,  j)  +  s{ij)  +  s{i,j  +  1)  +  s{i,j  -  1)]| 

This  formula  shows  that  the  discrete  Laplacian  operator  can  be  implemented  by 
subtracting  from  the  value  of  a  central  point  its  average  neighborhood.  Thus,  an  extended 
formula  [14]  can  be  written  as 

s{ij)  =  s{i,j)  +  k[s{i,j)-s{i,j)*h{i,j)],  (5) 

where  h{i,j)  is  a  discrete  averaging  filter,  and  *  denotes  convolution.  In  [14],  a 
equal-weighted  averaging  mask  was  used: 

r  1/N\  |x]<iV/2,|y|< 

1  0,  otherwise. 
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Another  way  to  extend  the  prototype  formula  [15]  came  from  the  idea  of 
Laplacian-of- Gaussian  hlter,  which  expands  Equation  (4)  into 

s{x,  y)  =  s{x,  y)  -  A:A  [s(a;,  y)  *  g{x,  y)]  =  s(x,  y)  -  k  [^(a;,  y)  *  Ag{x,  y)] ,  (6) 

where  g{x^y)  is  an  Gaussian  function,  and  Ag{x,y)  is  a  Laplacian-of-Gaussian  hlter. 

Finall}^,  we  mention  here  that  both  extensions  shown  in  Equations  (5)  and  (6)  are 
limited  to  a  single  scale. 

Inclusion  of  unsharp  masking  within  DDWT  framework 

We  shall  prove  that  unsharp  masking  with  a  Gaussian  lowpass  hlter  is  included  in  a 
dyadic  wavelet  framework  for  enhancement  by  considering  two  special  cases  of  linear 
enhancement. 

In  the  hrst  case,  transform  coefficients  of  channels  0  <  m  <  A  —  1  are  enhanced 
(multiplied)  by  the  same  gain  Go  >  1,  or  Gm  =Go>l,  0<m<A'  —  1.  The  system 
frequency  response  is  thus 

V{lj)  =  ^  GmCm{>^)  +  Cn{lo)  =  Go  ^  “  (^0  “  l)GAr(<^) 

m=0  m=0 

=  Go  —  (Go  —  l)GAf(c(;)  =  1  +  (Go  —  1)  [1  —  GA^(a;)]  . 

The  input-output  relationship  of  the  system  is  simply 

S[i]  =  s[i]  +  (Go  -  1)  {^[i]  -  *  c^r[^]}  .  (7) 

Since  is  approximately  a  Gaussian  lowpass  hlter.  Equation  (7)  may  be  seen  as  a 

1-D  counterpart  of  Equation  (5). 

In  the  second  case,  transform  coefficients  of  a  single  channel  p,  0  <  p  <  N  are  enhanced 
(multiplied)  by  a  gain  Gj,  >  1,  thus 

=  X]  +  GpCp{uj)  =  ^  Gm(<^)  +  (Gp  —  l)Gp(a;)  =  1  -T  (Gp  —  l)Gp(a;).  (8) 

m^p  171=0 

Using  the  hlter  class  of  Equation  (1),  the  input-output  relationship  of  the  system  (8) 
can  be  written  as 

s\i]  =  s[z]  -  (Gp  -  I)  ■  A  {^[z]  *  p[i]}  ,  (9) 

where  l3[i]  is  the  impulse  response  of  an  approximate  Gaussian  hlter.  Similarly,  Equation 
(9)  may  be  seen  as  a  I-D  counterpart  of  Equation  (6). 

The  inclusion  of  these  two  forms  of  unsharp  masking  clearly  demonstrates  the  hexibility 
and  versatility  of  a  dyadic  wavelet  framework. 
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2.1.4  Nonlinear  enhancement  by  functional  mapping 


Linear  enhancement  can  be  seen  as  a  mapping  of  wavelet  coefficients  by  a  linear  function 
Ern{x)  =  G-mX.  Therefore,  a  direct  extension  of  the  linear  enhancement  is  a  nonlinear 
mapping  function  Em{x). 

For  linear  enhancement,  selection  of  filters  G{(jj)  (and  thus  K{uij)  made  no  difference. 
However,  we  shall  see  that  the  selection  of  filters  is  critical  for  the  nonlinear  case. 


Nonlinear  enhancement  I:  Laplacian  filter. 

A  discrete  Laplacian  operator  can  be  implemented  by  the  filter 

l2 


G{lo)  =  -4 


sm 


(^)  ,  or  9[n]  =  {1,-2,!}, 


such  that  g[n]  *  s[n]  =  s[?r  +  1]  —  2s [n]  +  s[7r  —  1]. 

In  addition,  both  filters  iL(cu)  and  /L(cu)  can  also  be  symmetric, 


and 


K{^)  = 


H{u)  = 

G(u;) 


LO 

““ '  2/J 


T  2n 


2n-l  r 


=  E 

A  ^ 


1=0 


UJ 

cos  I- 


2/ 


Both  forward  and  inverse  filters  (0  <  m  <  iV  —  1)  can  be  derived  as 
F.ra{io)  =  —4  sin(2”^“^u;)  0m,2ra(<^)  =  — 4sin^  ^  4'”0m,2?i+2(^)  =  G{lXi)Era{oj)  (10) 


and 


1  2n— 1 

=  -0sn.2n(u7)-  E  H  (2™-la;)] 


2l 


r  f  ^  1 . 


1=0 


Note  that  the  forward  Liters  Em{oj)  (0  <  ?n  <  N)  can  be  interpreted  as  two  cascaded 
operations,  a  Gaussian  averaging  of  0m,2n+2(<^)  and  the  Laplacian  of  —4  |^sin(^)j  ,  while 
the  set  of  inverse  Liters  are  low-pass  Liters.  In  this  case,  wavelet  coefficients  may  be 

written  simply  as 

Wm[i]  =  A{s[f]  *  A„[z]} 

where  A  is  the  discrete  Laplacian  operator,  and  Am[f]  is  approximately  a  Gaussian  Liter. 
This  means  that  each  wavelet  coefficient  Wm[i\  is  dependent  on  the  local  contrast  of  the 
original  signal  at  i. 
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Figure  5:  (a)  E{x)  and  (b)  ^(x),  both  with  T  —  0.5  and  K  =  20. 


The  basic  guidelines  for  designing  a  nonlinear  enhancement  function  are: 


(1)  An  area  of  low  contrast  should  be  enhanced  more  than  an  area  of  high  contrast.  This 

is  equivalent  to  saying  that  small  values  of  Wm.\i]  should  have  larger  gains. 

(2)  A  sharp  edge  should  not  be  blurred. 

Such  functions  may  be  further  subjected  to  the  constraints  of: 

(1)  Monotonicity,  in  order  not  to  change  the  position  of  local  extrema,  nor  create  new 

extrema,  and 

(2)  Antisymmetry,  E{—x)  =  —E{x),  in  order  to  preserve  phase  polarity  for  “edge 

crispening”. 

A  simple  three-segment  function  that  satisfies  these  conditions  is  shown  in  Figure  5, 


£(,t)  = 


'  x~{K-l)T  ,]fx<-T 
Kx  ,  if  |a;|  <  r  ^  =  x  +  5(x) 

X  +  {K  —  l)T  ,  if  X  >  T 


(11) 


where  K  >  1  and 


8{x) 


'-{K-\)T,  ifx<-r, 
[K  —  l)a;,  if  |a:|  <  T, 

\k  -  i)r,  if  X  >  r. 


Thus,  an  enhanced  signal  can  be  written  as 

7V-1 

m=0 

N  N-1 

=  5[n]  =)=  fm[n]  *  im[n]  +  Y  (^N  *  ^mN)}  *  *mN 

m=0  m=0 
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Figure  6:  1-D  contrast  enhancement  by  four-level  dyadic  wavelet  analysis  with 
(a)  a  linear  operator  with  Kq  =  2.3,  and  (b)  a  nonlinear  operator  with  t  =  0.1  and  Ko  =  7- 

or, 

N-l 

s[?r]  =  5[n]  -  ^  {A  (^[n]  *  A^W)}  *  TmN-  (12) 

m=0 

We  point  out  that  the  formula  of  Equation  (12)  can  be  seen  as  a  multiscale  and 
nonlinear  extension  of  the  original  unsharp  masking  defined  by  Equation  (6). 

The  enhancement  operator  has  two  free  parameters:  threshold  Tm  and  gain  Km-  In 
our  exiDerimental  studies.  Km  =  Ko,  0  <  m  <  N  —  1,  and  Tm  =  t  x  max{|rnm[n]|},  where 
0  <  f  <  1  was  user  specified.  For  t  =  1.0,  wavelet  coefficients  at  levels  0  <  m  <  N  —  1  will 
be  multiplied  by  a  gain  of  Kq,  previously  shown  to  be  equivalent  to  unsharp  masking.  This 
means  our  nonlinear  algorithm  includes  unsharp  masking  as  a  subset.  Figure  6  shows  a 
numerical  example  comparing  linear  and  nonlinear  enhancement.  Note  the  lack  of 
enhancement  for  the  leftmost  edge,  for  the  case  of  the  linear  operator. 

We  argue  that  multiscale  unsharp  masking  as  defined  by  Equation  (12)  makes  an 
marked  improvement  over  traditional  techniques  in  two  respects: 

1.  The  fast  multiscale  (or  multimask)  decomposition  efficiently  searches  features  existing 
in  different  scales,  making  a  try-and-fail  strategy  of  window  selection  unnecessary. 

2.  The  nonlinear  algorithm  enhances  small  features  in  each  scale  without  blurring  edges 
of  larger  features,  making  simultaneous  enhancement  features  of  all  sizes  possible. 

Furthermore,  artifacts  possibly  created  by  a  nonlinear  enhancement  operator  can  be 
limited  by  careful  filter  selection  and  constraints.  For  example,  the  arguments  presented 
below  assure  that  no  new  extrema  will  be  created  within  each  channel. 


16 


1.  Filters  are  zero-phase.  No  spatial  shifting  exists  in  the  transform  space. 

2.  E{x)  is  a  monotonically  increasing  function  (Such  a  mapping  will  not  produce  new 
extrema  points.  Since  at  point  xq^  E  [/(xq)]  is  an  extrema  if  and  only  if  /(xq)  is  an 
extrema). 

3.  The  reconstruction  hlters  are  simply  zero-phase  smoothing  hlters. 


Nonlinear  enhancement  II:  gradient  filter. 

After  analyzing  the  Laplacian  filter,  a  natural  question  is,  “Are  there  any  other  possible 
filters  for  G{lo)T\  The  possible  candidates  are  limited  by  the  constraint 
G{uj)K{uj)  =  1  —  ||if((Xi)|p.  For  the  class  of  filters  defined  by  Equation  (1),  we  can  derive 
that 


G{lo)K{lo)  =  sin^ 


U) 


2n+p— 1  p 

E 

1=0 


U) 

cos  1  — 


1  2/ 


As  a  result,  the  only  other  meaningful  choice  of  filter  G{lo)  appears  to  be  a  gradient 
operator 

G{uj)  =  2ie"^2  sin  ,  or,  ^/[O]  =  l,g[l]  =  -1,  (13) 

such  that  g[i]  *  s[i]  =  ^[i]  —  —  1].  The  reasoning  comes  from  the  fact  that  derivatives 

higher  than  two  have  more  than  two  local  extrema  for  a  soft  edge.  Therefore,  nonlinear 
enhancement  of  derivatives  higher  than  two  may  produce  additional  edges  (artifacts). 

For  the  gradient  filter  G(a)),  we  selected  the  filter 


H[u)  = 


cos 


2n+l 


and  accordingly 


1 


2n 


^  ;=o 


LO 

cos  I  — 


1  2l 


We  then  derived  the  forward  filters 


E^,{lo)  =  G{u;)2^Qm,2n+2H  =  G(u;)A„,(cu) 


and  inverse  filters 


where 


r™(cu)  =  2™0 


m,2n+2 


21 
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is  a  low-pass  filter. 

In  this  case,  the  associated  wavelet  coefficients  may  be  written  as 

Wra[n]  =  V“  {^[n]  *  Xm[n]} 

where  is  a  discrete  backward  gradient  operator  characterized  by 
V“s[?7.]  =  s[n]  —  s[?7  —  1]. 

Finally,  the  enhanced  signal  may  be  written  as 

N-l 

('®W  *  ^mH)]  }  *  7m  w  +  5[n]  =1=  A7v[n]  *  7Ar[n]  (14) 

m=0 

where  V+  is  a  discrete  forward  gradient  operator  characterized  by  V+5[n]  =  s[n  +  1]  —  s[n] 
,  corresponding  to  e^^G{u)  in 

For  the  sake  of  discussing  artifacts,  let  ns  consider  the  channel  output  of  an  enhanced 
signal 

V+  jFlm  V“  (5m[n])  }  *  7mH 

where  5m  N  =  5[n]  *  A,„[?7]  and  7m  [?^]  is  a  zero-phase  low-pass  filter.  The  suspicious  part  is 

V+  {Em  V"  (sm[n])  }  , 
with  its  continuous  counterpart  being 

I  =  (x)  =  a[5;„(a;)]5"  (a;), 

where  Q!(x)  =  E'{x)  can  be  seen  as  a  gain  for  s"  (x).  Note  that  gain  a  is  based  on  5(„^(x) 
instead  of  s"  (x). 

The  proof  below  shows  that  only  a  linear  function  E{x)  =  Gx  (corresponding  to 
unsharp  masking)  can  guarantee  that  positions  of  extrema  of  a"  (x)  remain  unchanged. 

Proof.  Let  u{x)  =  a  [f'{x)]  /"(x),  then  u'{x)  =  a'(/')(/")^  +  «(/')/"'•  ^t  extrema 
point  .To  of  /"  ,  /'"(xo)  =  0,  thus  m'(xo)  =  a'  [/'(xq)]  [/"(xq)]^-  We  consider  a  particular 
signal  of  a  soft  edge,  /(x)  =  1/(1  -f  e“^^^),  for  which  f{x)  =  Pj  [2  cosh^(/9x)] , 
f"{x)  =  -/3^sinh(/3x)/cosh^(/?x)  and  /'"(x)  =  -(d^  [l  -  2sinh^(/?x)]  /  cosh^(/?x).  Local 
extrema  points  xq  satisfy  1  —  2sinh^(/?xo)  =  0,  f'{xo)  =  /9/3  and  [/”(xo)]^  =  4/3"‘/27.  For 
all  /3  0,  /'(xo)  0  and  [/"(xq)]^  0.  Therefore,  in  order  for  u'(xo)  =  0,  the  function 

a'{x)  must  satisfy  a'  [/'(xq)]  =  Q:'(/f/3)  =  0,V/3.  That  is  a[x)  =  constant,  or  E{x)  =  Gx. 

Change  in  position  of  extrema  points  may  result  in  undesirable  distortion  or  artifacts. 
Figure  7  shows  an  example  of  shifting  caused  by  a  hyperbolic  enhancement  function 


18 


Figure  7:  (a)  E  [5m(^)]  overlayed  with  5^(x)  (b)  a  [^^(a:)]  s'^{x)  overlayed  with  s'^{x). 

E{x)  =  A  ■  tanh(/Fx)  with  K  =  100.  For  other  nonlinear  functions,  new  local  extrema  may 
also  be  created. 

Intuitively,  the  problem  associated  with  the  gradient  filter  is  due  to  high-pass  filtering 
in  reconstruction.  The  significance  of  this  analytical  result  suggests  that  this  nonlinear 
approach  may  not  be  appropriate  for  other  wavelet  bases  if  their  reconstruction  (inverse) 
involves  high-pass  filtering. 

2.1.5  Noise-suppressed  enhancement 

A  structure  for  combined  denoising  and  enhancement 

The  nonlinear  enhancement  methods  proposed  previously  did  not  take  into  account  the 
presence  of  noise.  In  general,  noise  exists  in  a  digitized  image,  due  to  the  imaging  device 
and  quantization.  As  a  result  of  nonlinear  processing,  noise  may  be  amplified  and  may 
diminish  the  benefits  of  an  enhancement. 

Unfortunately,  denoising  is  a  very  difficult  problem  for  two  reasons.  Fundamentally, 
there  is  no  absolute  boundary  to  distinguish  a  feature  from  noise.  Even  if  there  are  known 
characteristics  of  a  certain  type  of  noise,  it  may  be  theoretically  impossible  to  completely 
separate  the  noise  from  features  of  interest.  Therefore,  most  denoising  methods  may  be 
seen  as  ways  to  suppress  very  high  frequency  and  incoherent  components  of  an  input  signal. 

A  very  simple  method  of  denoising  that  is  equivalent  to  low-pass  filtering  is  naturally 
included  in  our  dyadic  wavelet  framework.  That  is,  simply  discard  several  channels  of 
highest  resolution,  and  enhance  channels  confined  to  lower  frequency.  The  problem 
associated  with  this  linear  denoising  approach  is  that  edges  are  blurred  significantly,  and  is 
thus  not  suitable  within  a  contrast  enhancement  scheme.  Figure  10  (c)  shows  an  example 
of  this  approach.  In  order  to  achieve  edge-preserved  denoising,  more  sophisticated  methods 
based  on  a  framework  for  wavelet  analysis  were  proposed  in  the  literature.  Mallat  and 
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Hwang  [27]  connected  noise  behavior  to  singularities.  Their  algorithm  was  based  on  a 
multiscale  edge  representation.  The  algorithm  traced  modulus  wavelet  maxima  to  evaluate 
local  Lipschitz  exponents  and  deleted  maxima  points  with  a  negative  Lipschitz  exponent. 
Donoho  [25]  proposed  nonlinear  wavelet  shrinkage.  The  algorithm  shrinks  wavelet 
coefficients  towards  zero  based  on  a  level- dependent  threshold. 

The  wavelet  shrinkage  method  can  be  trivially  incorporated  into  our  nonlinear 
enhancement  framework  by  simply  adding  an  extra  segment  to  the  enhancement  function 
E{x)  of  Equation  (11). 

'  a,  _  (A^  _  l)r,  +  AT„  ,  if  .T  <  -Te 

K{x  +  r„)  ,  if  -  Te  <  X  <  -r„ 

E{x)  =  <  0  T  \x\  <  Tn  (15) 

K{x  -Tn)  ,ifTn<X<  Te 

a;  +  {K  -  1)7;  -  KT^  ,  ii  x>T, 

However,  there  are  two  arguments  which  favor  shrinking  gradient  coefficients  instead  of 
Laplacian  coefficients. 

First,  gradient  coefficients  exhibit  a  higher  signal  to  noise  ratio  (SNR).  For  any 
shrinkage  scheme  to  be  effective,  an  essential  property  is  that  the  magnitude  of  a  signal’s 
components  be  larger  than  that  of  existing  noise  (most  of  time).  It  is  thus  sensible  to 
define  the  SNR  as  the  maximum  magnitude  of  a  signal  over  the  maximum  magnitude  of 
noise.  For  example,  consider  a  soft  edge  model  f{x)  =  H/(l  +  H  >  0.  Its  first  and 

second  derivatives  are  f'{x)  =  Af^j  |^2 cosh^(/?a;)  and  f"{x)  =  —  sinh(/7x)/ cosh^(/3a;), 
with  magnitude  of  local  extrema  |/'(a:o)|  =  and  ]/"(a;o)|  =  2AI3'^ respectively. 

By  this  simple  model,  we  can  reasonably  assume  that  noise  is  characterized  by  relatively 
small  A  value  and  large  /5  value.  Clearly  that  gradient  coefficients  have  a  higher  SNR  than 
that  of  Laplacian  coefficients  because  of  less  contribution  of  /?.  Figures  8  (b)  and  (c)  show 
first  and  second  derivatives,  respectively,  for  an  input  signal  (a)  with  two  distinct  edges. 

Second,  boundary  contrast  is  not  affected  by  shrinking  gradient  coefficients.  As  shown 
in  Figure  8,  coefficients  aligned  to  the  boundary  of  edges  are  the  local  extrema  in  the  case 
of  the  hrst  derivative  (gradient),  and  zero  crossings  in  the  case  of  the  second  derivative 
(Laplacian).  For  a  simple  point-wise  shrinking  operator,  there  is  no  way  to  distinguish  the 
B  points  from  the  A  points.  As  a  result,  both  regions  around  the  B’s  and  A’s  are  shrunken, 
and  the  discontinuity  in  B  will  sacrifice  boundary  contrast. 

In  the  previous  section,  we  argued  that  nonlinear  enhancement  should  be  performed  on 
Laplacian  coefficients.  Therefore,  in  order  to  incorporate  denoising  into  our  enhancement 
algorithm,  we  need  to  split  the  Laplacian  operator  into  two  cascaded  gradient  operators. 
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Figure  8:  (a)  Signal  with  two  edges,  (b)  1st  derivative  (gradient), 
(c)  2nd  derivative  (Laplacian).  (d)  Shrunken  2nd  derivative. 


Figure  9:  Incorporating  wavelet  shrinkage  into  an  enhancement  framework  (level  one  shown). 


Note  that 


=  —4  sin  ^2™  = 


,  if  m  =  0, 

6'(i(2™“^c<;)]^  ,  otherwise. 


(16) 


where  Gdiyj)  =  2j  sin(i<;). 

Denoising  by  wavelet  shrinkage  [25]  can  then  be  incorporated  into  this  structure  as 
illustrated  in  Figure  9,  where  the  shrinking  operator  can  be  written  as 


C{x)  =  sign{x)  ■ 


x\  -Tn  ,if  |a;l  >  r„, 
0  ,  otherwise. 


Note  that  the  shrinking  operator  is  a  piece-wise  linear  and  monotonically  non-decreasing 
function.  In  practice,  this  will  not  introduce  artifacts. 

Finally,  a  denoised  and  enhanced  signal  can  be  written  overall  as 


sm 


7V-1 

-  X]  {v+  [Cm  (V“  (^[n]  *  A„[n]))]|  *  7^[n]  -k  s[n]  *  *  7iv[n]. 

m=0 


21 


Threshold  estimation  for  denoising 

The  threshold  Tn  is  a  critical  parameter  in  the  shrinking  operation.  For  a  white  noise 
model  and  orthogonal  wavelet,  Donoho  [25]  suggested  a  formula  of  T„  =  ^2  \og{N)a  j ^/N ^ 
where  N  is  the  length  of  a  input  signal  and  a  is  the  standard  deviation  of  wavelet 
coefficients.  However,  the  dyadic  wavelet  we  applied  is  not  an  orthogonal  wavelet. 
Moreover,  in  our  2-D  applications,  a  shrinking  operation  is  applied  to  magnitudes  of 
gradient  coefficients  instead  of  wavelet  coefficients  themselves.  Therefore,  the  threshold 
estimation  method  proposed  in  [28]  for  edge  detection  may  be  more  suitable. 

In  our  shrinking  operation,  the  magnitudes  of  the  gradient  of  a  Gaussian  low-passed 
signal  are  modified.  As  pointed  out  in  [28],  for  white  Gaussian  noise,  the  probability 
distribution  function  of  the  magnitudes  of  gradient  is  characterized  by  the  Rayleigh 
distribution: 

I"  r|e-(mAF/2  ^>0, 

Pr\\/xf\\{m)  =  < 

0  ,m  <  0. 

To  estimate  t/,  a  histogram  (probability)  of  ||A/||  was  computed,  and  then  iterative  curve 
fitting  was  applied.  Under  this  model,  the  probability  p  of  noise  removal  for  a  particular 
threshold  r  can  be  calculated  [28]  by 

^  /o~Pr||A/||(m)dm’ 

and  thus  r  =  21n(l  —  p)  p.  For  p  =  0.999,  r  =  3.7//. 

Figure  10  compares  different  existing  approaches.  In  (b),  we  observed  that 
enhancement  without  any  denoising  results  in  annoying  background  noise.  In  (c),  edges 
were  smeared  and  broadened  by  low-pass-enhancement  combination.  Only  in  (d),  with 
wavelet  shrinkage,  we  were  able  to  achieve  the  remarkable  result  of  denoising  and  contrast 
enhancement  simultaneously. 

To  demonstrate  the  denoising  process.  Figure  11  (a)  and  (b)  shows  both  nonlinear 
enhancement  of  wavelet  coefficients  without  and  with  denoising,  respectively,  for  the  input 
signal  shown  in  Figure  10  (a).  Figure  11  (c)  shows  the  curve-fitting  for  threshold 
estimation. 

2.1.6  Two  dimensional  extension 

For  image  processing  applications,  the  one  dimensional  structures  discussed  previously 
were  simply  extended  to  two  dimensions.  We  adopted  the  method  proposed  by  Mallat  [26], 
shown  in  Figure  12,  where  filter  L{lo)  =  lilRhdL^  K{^)  ^iicl  G{io)  are  the  same 

filters  used  in  the  1-D  case. 
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(a) 


(d) 


Figure  10:  (a)  Noisy  input  signal  (contaminated  by  white  Gaussian  noise). 

(b)  Nonlinear  enhancement  without  denoising,  Gm  =  10,  =  4,  t  =  0.1. 

(c)  Nonlinear  enhancement  of  levels  2-3,  Gm  =  10,  t  =  0.1;  levels  0-1  zeroed  out; 
(d)  Nonlinear  enhancement  with  adaptive  wavelet  shrinkage  denoising,  Gm  =  10,  =  4,  t 


Figure  11:  Column  (a),  Enhanced  wavelet  coefficients  without  denoising. 
Column  (b),  Enhanced  wavelet  coefficients  with  adaptive  thresholding  Tn  —  d.h??. 
Column  (c),  The  magnitude  distribution  and  curve-fitting. 

(Row  1  through  4  corresponds  to  levels  1  to  4.) 


Figure  12:  Two  dimensional  dyadic  wavelet  transform  (two  levels  shown). 
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•Gd(cOx); 


Figure  13:  Denoising  and  enhancement  for  the  2-D  case  (level  one  shown). 

However,  experimental  we  observed  that  if  we  simply  modify  the  two  oriented  wavelet 
coefficients  independently,  we  introduced  orientation  distortions.  One  way  to  avoid  this 
disastrous  artifact  is  to  apply  a  denoising  operation  on  the  magnitude  of  gradient 
coefficients,  and  apply  a  nonlinear  enhancement  operation  on  the  sum  of  the  Laplacian 
coefficients,  as  shown  below  in  Figure  13.  For  the  two  oriented  gradient  coefficients  wxi 
and  wyi ,  the  magnitude  M  and  phase  P  were  computed  as  M  =  ^Jwxl  +  wyj  and 
P  =  arctan(t(;yi/roxi),  respectively.  The  denoising  operation  was  then  applied  to  M, 
obtaining  M' .  The  denoised  coefficients  were  then  simply  restored  as  wx[  =  M'  *  cos(P) 
and  wy[  =  M'  *  sin(P),  respectively.  For  the  enhancement  operation,  notice  that  the  sum 
of  two  Laplacian  components  is  isotropic.  Therefore,  we  may  compute  the  sum  of  the  two 
Laplacian  components  S  =  wx^  +  wy2  and  C  =  WX2I S.  A  nonlinear  enhancement  operator 
is  then  applied  to  only  S',  producing  S' .  Thus,  the  two  restored  components  are 
wx'2  =  S'  *  C  and  wy'2  =  S''  *  (1  —  C). 


2.1.7  Experimental  results  and  comparisons 

In  this  section,  we  present  sample  experimental  results  and  compare  them  with  existing 
state-of-the-art  techniques. 

Figure  14  (a)  shows  a  synthetic  image  with  three  circular  “bumps”  and  added  white 
noise.  The  enhancement  results  shown  in  (b)  and  (c)  demonstrate  unwanted  noise 
amplification.  Moreover,  note  that  histogram  equalization  processing  alters  the  object’s 
boundary.  However,  the  result  shown  in  (d)  accomplished  by  dyadic  wavelet  analysis 
produces  a  clearer  image  without  orientation  distortion. 

Figure  15  (a)  shows  an  original  dense  mammogram  image  with  a  central  mass.  The 
boundary  of  the  mass  in  the  the  enhanced  image  was  well  defined  and  made  clear 
penetration  of  spicules  into  the  mass. 

To  study  the  efficacy  of  our  algorithm,  we  blended  mathematical  phantom  features  into 
clinically  proved  cancer  free  mammograms.  Figures  16  (a)  and  (b)  show  mathematical 
phantom  features  blended  into  M48  and  M56  (resulting  in  Figure  17  (a)  and  Figure  18 
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(c)  (d) 

Figure  14:  (a)  Noisy  image  (white  Gaussian  noise  contaminated). 

(b)  Histogram  equalized,  (c)  Nonlinear  enhancement  by  Beghdadi  and  Negrate’s  algorithm, 
(d)  Nonlinear  enhancement  with  adaptive  wavelet  shrinkage  denoising,  Gm  =  20,  =  4,  t  =  0.1. 
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0)  (b) 


Figure  15:  (a)  Oringinal  mammogram  image  M73.  (b)  Xonlinear  enhancement  with  adaptive 

wa^■elet  shrinkage  denoising.  G’,,-,  =  20.  .V  =  5.  /  =  0.1. 


Figure  16:  (a)  Phantom  features  blended  into  MdS.  (b)  Phantom  features  blended  into  M56. 
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(a)),  respectively. 

Figure  17  (a)  shows  a  dense  mammogram  with  blended  phantom  features,  and  (b) 
shows  an  image  processed  by  our  nonlinear  method.  The  enhanced  image  delineates  well 
the  boundary  (uncompressed  areas)  of  the  breast  and  its  structure.  The  phantom  features 
were  also  well  enhanced.  Figure  18  (a)  shows  a  dense  mammogram  with  blended  phantom 
features,  and  (b)  shows  the  image  enhanced. 

2.1.8  Summary 

We  established  connections  between  dyadic  wavelet  enhancement  algorithms  and 
traditional  unsharp  masking.  We  proved  that  two  cases  of  linear  enhancement  were 
mathematically  equivalent  to  traditional  nnsharp  masking  with  Gaussian  low-pass  filtering. 
We  designed  a  methodology  for  nonlinear  enhancement  with  a  simple  nonlinear  function  to 
overcome  the  dynamic  range  requirement  nsually  associated  contrast  enhancement  of 
digital  radiographs.  By  careful  selection  of  wavelet  filters  and  enhancement  function,  we 
showed  that  artifacts  can  be  eliminated.  An  additional  advantage  of  our  simple 
enhancement  function  is  that  it  includes  traditional  unsharp  masking  as  a  subset.  We 
showed  how  an  edge-preserved  denoising  stage  (wavelet  shrinkage)  can  be  appropriately 
incorporated  into  our  contrast  enhancement  framework,  and  introduced  a  method  for 
adaptive  threshold  estimation.  We  showed  how  denoising  and  enhancement  operations 
should  be  carried  out  for  two  dimensional  images  to  avoid  orientation  distortions. 

Our  future  research  plan  shall  include  the  systematic  study  of  gain  and  threshold 
parameters  for  the  nonlinear  enhancement.  In  addition,  in  the  next  year  we  shall  seek 
localized  and  complex  nonlinear  methods  to  improve  the  performance  of  our  existing 
algorithm. 

2.2  XMam:  An  Image  Editor  for  Mammography 

2.2.1  Introduction 

Recent  estimates  show  that  breast  cancer  is  the  most  frequently  diagnosed  malignancy  in 
the  United  States,  accounting  among  women  for  32%  of  all  cancers  detected  and  18%  of  all 
cancer  deaths  [3].  Although  mammography  is  widely  recognized  as  the  best  method  for  the 
detection  of  breast  cancer,  10%-30%  of  women  who  have  breast  cancer  and  undergo 
mammography  have  negative  mammograms,  and  in  approximately  two-thirds  of  these 
false-negative  mammograms  the  radiologist  failed  to  identify  cancers  that  were  visible  upon 
retrospective  review  of  the  radiographs  [29].  In  [30],  Vyborny  suggests  that  part  of  the 
problem  resides  in  human  search  performance,  a  conclusion  he  bases  on  the  following 
findings:  (1)  Lack  of  prior  knowledge  of  the  location  of  an  abnormality  decreases  the 
likelihood  that  it  will  be  detected,  (2)  In  chest  radiography,  a  large  portion  of  the  image  is 
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not  sampled  by  foveal  vision  during  the  radiologist’s  visual  search,  an  effect  that  is  also 
observed  in  mammography,  although  to  a  lesser  extent  because  of  the  smaller  size  of  the 
images,  and  (3)  In  searching  for  lung  nodules,  observers  exhibit  longer  visual  dwell  times 
on  a  considerable  portion  of  missed  radiologic  abnormalities  than  on  normal  areas,  a 
phenomenon  that  is  significant  enough  to  allow  for  computer-based  visual  aids  to  reduce 
false-negative  detections.  Both  [29]  and  [30]  point  out  that  replicated  mammography 
readings  (by  two  radiologists  or  a  radiologists  and  a  technologists)  may  increase  selectivity 
and  decrease  the  number  of  false-negative  mammograms.  One  of  the  goals  of 
computer-aided  diagnostic  schemes  is  to  aid  the  radiologist  in  the  search  for  lesions  by 
indicating  locations  of  suspicious  abnormalities  in  mammograms. 

When  a  suspicious  region  is  detected,  the  radiologist  must  visually  extract  features 
from  the  finding  and  then  either  decide  whether  the  abnormality  is  malignant  or  benign,  or 
recommend  additional  screening,  follow-up  or  biopsy.  Although  there  are  general  guidelines 
to  differentiate  benign  from  malignant  breast  lesions,  a  considerable  number  of  lesions  are 
misclassified.  For  example,  more  than  half  of  false-negative  mammograms  result  from 
decisions  by  radiologists  that  the  finding  is  either  normal  or  benign  [30],  while  only 
10%-20%  of  masses  in  patients  referred  for  surgical  breast  biopsy  are  malignant  [29].  In 
[30],  Vyborny  suggests  that  part  of  the  problem  resides  in  the  merging  of  multiple  image 
features  to  characterize  abnormalities.  Although  general  radiologist  can  accurately  and 
reproducibly  extract  individual  image  features  of  importance  in  mammography,  they  are 
usually  outperformed  by  experienced  mammographers  in  combining  multiple  features  into 
a  correct  diagnosis  [30].  One  objective  of  computer-aided  diagnostic  schemes  is  to  extract 
and  analyze  the  characteristics  of  lesions  to  aid  the  radiologist  in  their  classification. 

Fundamental  to  the  success  of  a  computer-aided  diagnostic  methodology  is  the 
development  of  a  large  database  of  cases  diagnosed  as  having  signs  of  cancer  based  on 
mammographic  screening,  and  later  conhrmed  or  denied  as  positive  by  biopsy.  This 
database  will  serve  as  the  ground  truth  for  the  development  and  verification  of  our 
computer  algorithms  for  the  detection  of  breast  cancer.  In  this  section  of  our  report  we 
describe  the  first  version  of  a  computer  image  editor  that  provides  the  necessary 
functionality  to  compile  such  database  of  cases. 

2.2.2  Functionality 

XMam  is  a  software  tool  that  allows  radiologist  to  interactively  indicate  on  digitized 
mammograms  regions  diagnosed  as  having  signs  of  cancer.  XMam’s  user  interface  has  been 
designed  to  fit  the  needs  of  mammographic  screening  and  consists  of  four  image  panels  to 
allow  for  the  simultaneous  display  of  the  craniocaudal  (CC)  and  mediolateral  oblique 
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(MLO)  views  of  the  left  and  right  breasts  (the  four  accepted  views  for  breast  cancer 
screening  in  the  United  States).  As  shown  in  Figure  19,  XMam’s  interface  is  made  up  of 
two  top  windows  which  logically  group  the  images  corresponding  to  the  CC  and  MLO 
views.  These  windows  consist  of  left  and  right  image  panels  equipped  with  vertical  and 
horizontal  scrollbars.  The  CC  window  is  the  application’s  main  window  and  provides,  in 
addition  to  its  two  image  panels,  three  pull-down  menus  and  one  radio-box.  In  overview, 
XMam’s  functionality  is  supported  by  three  logical  modules;  (1)  an  image  module,  (2)  a 
data  module  and  (3)  a  drawing  module.  In  the  next  three  sections  we  describe  each  of 
these  modules  in  detail. 

Image  module 

The  image  module  manages  image  retrieval  and  display  operations.  This  module 
provides  a  user  with  the  necessary  functionality  for  opening  a  patient’s  image  file  and 
displaying  its  contents  in  the  corresponding  windows.  In  the  context  of  XMam,  an  image 
file  refers  to  the  logical  grouping  of  the  four  standard  radiographs  obtained  from  a  single 
mammographic  screening.  If  for  any  reason,  one  or  more  images  are  missing  from  the 
patient’s  image  file,  the  image  module  will  leave  blank  the  corresponding  image  panel(s). 

A  user  may  access  the  image  module  through  the  Open  and  Close  buttons  located  in 
the  Image  pull-down  menu  on  XMam’s  main  window.  The  layout  of  this  menu  as  well  as 
other  pull-down  menus  is  shown  in  Figure  20.  A  point-and-click  operation  on  the  Open 
button  of  the  Image  pull-down  menu  pops  up  a  File  Image  Dialog  as  the  one  shown  in 
Figure  21,  and  allows  a  user  to  retrieve  a  patient’s  image  file  from  the  group  of  patients 
listed  in  its  directory.  A  user  may  close  a  patient’s  image  file  either  by  opening  a  new  file 
from  the  File  Image  Dialog  or  with  a  point-and-click  operation  on  the  Close  button  of 
the  Image  pull-down  menu. 

Data  module 

This  module  manages  data  update  and  retrieval  operations.  It  provides  a  user  with  the 
necessary  functionality  for  loading  and  saving  a  complex  data  structure  that  indicates 
which  image  regions  were  registered  by  a  radiologist  as  having  signs  of  cancer.  A  complete 
description  of  this  data  structure  is  described  in  detail  below. 

A  user  may  access  the  data  module  through  the  Load  and  Save  buttons  located  in  the 
Data  pull-down  menu  on  XMam’s  main  window.  A  point-and-click  operation  on  the  Load 
button  of  the  Data  pull-down  menu  pops  up  a  Load  Dialog,  as  shown  in  Figure  22(a),  and 
allows  a  user  to  retrieve  a  patient’s  data  structure.  Similarly,  a  point-and-click  operation 
on  the  Save  button  of  the  Data  pull-down  menu  pops  up  a  Save  Dialog  as  shown  in 
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Figure  19:  XMam’s  user  interface:  (a)  CC  window,  (b)  MLO  window. 
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Figure  20:  XMam’s  pull-down  menus:  (a)  Image  menu,  (b)  Data  menu 
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Figure  22:  XMam’s  File  Data  Dialogs:  (a)  Load  Dialog,  (b)  Save  Dialog. 

Figure  22(b),  and  allows  a  user  to  store  the  data  structure  that  contains  the  specific 
regions  indicated  by  a  radiologist  as  having  signs  of  cancer. 

Drawing  module 

The  drawing  module  manages  the  drawing  operations  of  XMam.  This  module  provides 
a  user  with  the  necessary  functionality  for  outlining  image  regions.  This  is  accomplished 
with  the  aid  of  a  computer  mouse  by  enclosing  regions  of  interest  with  polygons. 

Before  describing  in  detail  XMam’s  drawing  operations,  we  introduce  some  terminology 
that  will  simplify  our  discussion.  In  XMam,  an  active  polygon  group  refers  to  a  group  of 
polygons  on  which  a  user  is  performing  drawing  operations.  At  any  point  in  time  there  can 
be  at  most  one  such  group.  In  general,  a  polygon  group  becomes  active  when  a  user 
performs  a  point-and-click  operation  on  any  of  its  polygons.  A  polygon  that  is  being  drawn 
is  considered  as  active  and  forms  a  polygon  group  on  its  own.  Polygon  groups  with  more 
than  one  polygon  are  constructed  by  selecting  several  polygons  and  grouping  them  together 
by  a  group  operation.  A  polygon  group  can  be  viewed  conceptually  as  an  unordered  tuple 
of  elements.  Tuple  elements  may  be  single  polygons  or  tuples  themselves.  This  tuple 
feature  allows  a  user  to  impose  a  complex  logical  structure  into  a  polygon  group.  As  an 
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Figure  23:  Outlining  of  a  region  of  interest  on  XMam. 


example,  suppose  that  a  and  b  are  two  polygons  enclosing  regions  of  interests  corresponding 
to  a  single  lesion  visible  on  a  given  mammographic  view,  say  the  left  craniocaudal  (LCC). 
Then,  a  and  b  can  be  grouped  together  on  a  single  tuple  (a,  b)  to  indicate  that  they  are 
logically  related.  Also,  suppose  that  c  is  a  polygon  enclosing  a  region  related  to  the  same 
lesion  but  visible  on  the  complementary  left  mediolateral  view  (LMLO).  To  complete  the 
lesion’s  description,  polygon  c  should  be  grouped  with  tuple  (a,  b)  to  yield  a  complex  group 
((a,  5),  c).  An  important  advantage  of  this  tuple  structure  is  that  image  regions  that  belong 
to  the  same  lesion  but  that  are  visible  on  different  views  can  be  logically  grouped  together. 
Notice  that  this  information  may  be  useful  for  algorithms  that  rely  on  the  simultaneous 
fusion  of  all  four  views  for  the  detection  of  breast  cancer.  Another  advantage  of  having 
polygon  groups  is  that  when  a  user  selects  any  polygon,  all  other  polygons  that  form  part 
of  its  group  become  active.  This  provides  immediate  feedback  as  to  which  image  regions 
are  logically  related.  In  addition  this  functionality  can  later  be  used  for  training  new 
mammographers  and  technicians  with  previously  diagnosed  cases  of  breast  cancer. 

There  are  several  drawing  operations  a  user  may  engage  at  any  given  point  in  time. 

The  simplest  operation  is  that  of  drawing  a  polygon.  Figure  23  shows  several  of  the  steps 
followed  by  a  user  to  enclose  a  suspicious  mass  on  a  digitized  mammogram.  Recall  that  a 
polygon  that  is  being  drawn  is  considered  an  active  polygon  group,  a  state  that  is  denoted 
in  XMam  by  highlighting  the  vertices  of  the  active  polygons  with  square  anchors.  Figure 
23  also  shows  that  there  are  two  classes  of  anchors  in  XMam,  active  and  inactive.  An 
active  anchor  is  displayed  as  a  solid  square  and  denotes  a  polygon  vertex  on  which  the  user 
may  perform  drawing  operations.  On  the  other  hand,  an  inactive  anchor  is  displayed  as  a 
clear  square  and  cannot  be  operated  on  until  it  is  in  the  active  state.  At  any  point  in  time 
there  can  be  at  most  one  active  anchor  within  a  polygon  group.  An  anchor  becomes  active 
when  a  user  performs  a  point-and-click  operation  on  it.  When  an  anchor  is  active,  a  user 
may  adjust  the  anchor  position,  delete  the  anchor,  or  adjust  the  polygon  position. 

XMam  provides  four  different  mouse  cursors  to  indicate  its  current  state.  We  will  refer 
to  these  cursors  as  crosshair,  arrow,  dotbox,  and  fleur.  A  crosshair  cursor  indicates  that  a 
user  may  start  drawing  a  polygon  at  a  given  mouse  position.  If  a  user  places  the  mouse 
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cursor  on  top  of  any  line  segment  of  an  inactive  polygon,  the  cursor  will  change  from  a 
drawing  cursor  to  an  arrow  cursor  to  indicate  that  the  user  may  select  the  inactive  polygon 
by  a  point-and-click  operation.  Similarly,  if  a  user  places  the  mouse  cursor  on  top  of  any 
line  segment  of  an  active  polygon,  the  cursor  will  change  from  a  drawing  cursor  to  a  dotbox 
cursor  to  indicate  that  the  user  may  add  an  additional  anchor  to  the  polygon  at  the 
current  mouse  position.  Also,  if  a  user  places  the  mouse  cursor  on  top  of  any  vertex  of  an 
active  polygon,  the  mouse  cursor  will  change  from  a  drawing  cursor  to  a  fleur  cursor  to 
indicate  that  the  user  may  adjust  the  anchor  or  polygon  position  by  a  point-and-click 
operation.  The  operation  to  be  performed  depends  on  which  mouse  button  is  pressed 
during  the  point-and-click  operation.  Adjustment  of  the  anchor  or  polygon  position  may  be 
carried  out  until  the  user  releases  the  mouse  button.  In  addition,  a  user  may  delete  an 
active  anchor  by  depressing  the  delete  key. 

Recall  that  there  can  be  at  most  one  active  polygon  group  at  any  given  time.  We  also 
mentioned  that  a  polygon  group  with  more  than  one  polygon  can  be  constructed  by 
selecting  several  polygons  and  performing  a  group  operation.  Selecting  a  polygon  makes  it 
active,  therefore  we  needed  a  special  operation  to  override  the  constraint  that  at  most  one 
polygon  group  may  be  active  at  any  point  in  time.  This  is  accomplished  in  XMam  by 
depressing  the  shift  key.  Polygon  groups  selected  while  depressing  the  shift  key  become 
active  and  can  be  grouped  together  by  a  group  operation.  If  the  group  operation  is  not 
performed  and  the  shift  key  is  released,  then  selecting  a  new  polygon  will  make  it  active 
and  all  other  polygons  inactive. 

XMam  not  only  allows  a  user  to  select  regions  of  interests  but  also  provides  the 
functionality  to  classify  them  under  five  distinct  lesion  classes:  (1)  well-defined  masses,  (2) 
ill-defined  masses,  (3)  spiculated  lesions,  (4)  architectural  distortions  and  (5) 
microcalcifications.  Figure  24  shows  a  radio-box  containing  five  radio-buttons,  each 
associated  with  a  distinct  lesion  class.  Through  this  menu,  a  user  may  define  a  class  under 
which  any  polygons  drawn  on  the  image  panels  are  to  be  classified  and  displayed.  In  other 
words,  the  lesion  menu  provides  a  window  into  the  lesion  space  and  allows  the  user  to 
display  only  those  polygons  that  are  related  to  a  given  lesion  type.  In  the  current  version 
of  XMam,  a  user  is  only  allowed  to  display  at  any  given  time  polygons  belonging  to  a  single 
lesion  type.  This  limitation  was  imposed  by  our  existing  Sun  workstation  with  8-bit 
frame-buffers.  Although  there  exist  several  software  solutions  for  the  problem,  our  solution 
is  a  hardware  upgrade  to  a  24-bit  frame-buffer.  A  temporary  software  solution  was  not 
implemented  because  we  anticipate  to  port  XMam  to  a  24-bit  frame-buffer  platform  in  the 
near  future  and  determined  that  clouting  (patching)  the  software  with  a  partial  solution 
would  be  detrimental  to  the  real-time  performance  of  the  tool. 
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Figure  24:  XMam’s  lesion  menu. 


2.2.3  Data  structure 

This  section  describes  the  data  structures  used  by  XMam  to  support  drawing  and  grouping 
operations.  As  discussed  above,  a  polygon  group  can  be  viewed  as  an  unordered  tuple  of 
elements,  where  tuple  elements  may  be  single  polygons  or  tuples  themselves.  In  order  to 
support  this  recursive  structure,  XMam  represents  a  tuple  by  a  linked  list  of  cells,  where 
cells  may  point  to  single  polygons  or  other  linked  lists.  Cells  that  point  to  single  polygons 
are  referee!  to  as  polygon  cells,  whereas  cells  that  point  to  other  linked  lists  are  refered  to 
as  group  cells.  Figure  25  shows  an  instance  of  each  of  these  two  cell  types  with  their 
corresponding  data  fields  (the  field  used  to  identify  the  cell  type  was  omitted  for 
simplicity.)  For  any  given  cell  c,  PPrevCell  and  PMextCell  are  pointers  to  the  cells 
preceding  and  succeeding  c.  For  a  polygon  cell  p,  PArray  serves  as  a  pointer  to  an  array 
that  sequentially  lists  the  X  and  Y  coordinates  of  the  polygon  associated  with  p,  Panel 
indicates  the  panel  in  which  p’s  polygon  will  be  displayed  (LCC,RCC,LMLO,RMLO),  TOL 
indicates  the  type  of  lesion  p’s  polygon  belongs  to  (well-defined  mass,  ill-defined  mass, 
etc.),  NPoints  indicates  the  number  of  vertices  of  p’s  polygon,  and  FPoint  indicates  the 
index  (if  any)  of  the  active  anchor  in  p’s  polygon.  For  a  group  cell  g,  PGFirst  and  PGLast 
are  pointers  to  the  first  and  last  cells,  respectively,  of  the  linked  list  associated  with  p’s 
group.  By  including  the  Panel  and  TOL  fields  on  each  polygon  cell  we  may  create  polygon 
groups  consisting  of  polygons  from  different  panels  and  having  different  lesion  types. 
Although  our  current  user  interface  only  allows  a  user  to  form  groups  within  a  specific 
lesion  type,  our  data  structure  already  supports  grouping  across  different  lesion  types.  This 
additional  grouping  feature  is  provided  in  view  of  our  plan  to  port  XMam  to  a  24-bit 
framebuffer  platform  and  use  colors  to  indicate  the  lesion  tyjae  of  a  polygon. 

To  keep  track  of  all  polygon  groups,  XMam  compiles  them  into  a  single  linked  list.  This 
list  is  accessed  through  a  handle  consisting  of  two  fields,  PFirst  and  PLast,  that  serve  as 
pointers  to  the  first  and  last  cells  of  the  list,  respectively.  In  Figure  26  we  show  how  the  list 
handle,  the  polygon  cells  and  the  group  cells  are  used  together  to  represent  the  polygon 
group  ((a,  6),  c)  described  in  the  previous  section.  We  will  refer  to  this  complex  data 
structure  as  the  region  structure.  A  more  complex  example  consisting  of  two  polygon 
groups,  (a,  6,  c)  and  ((d,  e),  (/,(/)),  is  shown  in  Figure  27. 

The  region  structure  can  be  easily  traversed  by  following  the  PNextCell  pointers  of  its 
linked  list.  Every  time  a  group  cell  is  found,  it  is  recursively  expanded  and  its  underlying 
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list  traversed.  An  example  of  this  operation  occurs  when  a  user  changes  the  type  of  lesion 
to  be  displayed  and  XMam  traverses  its  region  structure,  clearing  the  polygons  belonging 
to  the  previous  lesion  type  and  drawing  the  polygons  corresponding  to  the  new  lesion  type. 

In  addition  to  the  region  structure,  XMam  uses  an  auxiliary  structure  called  the  group 
structure  to  keep  track  of  the  active  polygon  group.  Figure  28  shows  an  instance  of  the 
group  structure  consisting  of  a  list  handle  and  four  list  cells.  The  LFirst  and  LLast  fields 
serve  as  pointers  to  the  first  and  last  cells  of  the  list,  respectively.  Each  list  cell  in  turn 
points  to  a  polygon  group  selected  by  a  user.  This  list  usually  consists  of  a  single  entry 
pointing  to  the  active  polygon  group  but  may  contain  several  entries  when  a  user  selects 
several  polygon  groups  for  a  grouping  operation.  The  other  three  fields  in  the  list  handle 
function  as  different  levels  of  addressing  for  the  polygon  that  contains  the  active  anchor. 
The  LCurr  field  points  to  the  list  cell  that  points  to  the  polygon  group  that  contains  the 
polygon  cell  with  the  active  anchor.  The  PGCell  field  points  to  the  group  that  directly 
contains  the  polygon  cell  with  the  active  anchor.  And  the  PPCell  field  directly  points  to 
the  polygon  cell  that  contains  the  active  anchor.  The  reason  for  having  these  three  distinct 
pointers  is  justified  below. 

We  clearly  require  a  pointer  to  the  polygon  cell  with  the  active  anchor  to  support 
possible  edit,  insert  and  delete  operations.  In  addition,  if  a  user  deletes  a  polygon  from  a 
polygon  group  consisting  of  two  polygons,  then  the  group  cell  becomes  unnecessary  and 
needs  to  be  deleted.  The  PGCell  field  provides  efficient  access  to  this  cell.  Similarly,  if  the 
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Figure  26:  Structure  for  polygon  group  ((a,  6),  c). 
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Figure  27:  Structure  for  polygon  groups  (a,  6,c)  and  ((d,  e),  {f,g))- 


Figure  28:  Group  structure. 


polygon  group  containing  the  active  anchor  consists  of  a  single  polygon  and  a  user  deletes 
this  polygon,  then  the  list  cell  becomes  nnnecessary  and  needs  to  be  deleted.  The  LCurr 
held  provides  efficient  access  to  this  cell. 

Although  XMam  uses  several  other  data  structures,  the  region  and  group  structures 
described  above  comprise  the  major  mechanisms  used  by  XMam  for  managing  drawing  and 
grouping  operations.  Even  though  these  data  structures  fully  support  the  functionality 
described  in  the  previous  section,  we  anticipate  that  future  versions  of  XMam  shall  yield 
cleaner  and  more  efficient  data  structures. 

2.2.4  Hardware 

In  this  section  we  briefly  describe  our  current  hardware  platform.  The  development 
environment  for  XMam  is  a  dual-screen  SPARCstation  LX  with  8-bit  framebuffers,  64 
Mbytes  of  memory  and  320  MBytes  of  swap  space.  The  platform  currently  runs  Sun  OS 

4.1.3  and  XWindows  X11/R5,  and  includes  the  recently  released  Xvan  server.  This  new 
server  allows  us  to  treat  both  screens  as  a  single  logical  device  and  to  freely  move  (drag) 
windows  from  one  screen  to  the  other. 

2.2.5  Summary 

We  have  described  above  the  functionality  and  data  structures  of  a  computer  image  editor 
that  can  be  used  by  radiologists  to  indicate  on  digitized  mammograms  regions  diagnosed  as 
probable  cancer.  This  tool  shall  facilitate  the  compilation  of  a  large  database  of  diagnosed 
cases  to  be  used  as  ground  truth  for  the  development  and  verification  of  our  computer 
algorithms  for  the  detection  of  breast  cancer. 

In  the  next  year,  XMam  will  progresively  evolve  into  a  state-of-the-art  mammography 
workstation  as  new  software  is  developed  and  hardware  improvements  are  incorporated 
into  our  system. 

2.3  Local  Feature  Analysis  via  Interval  Wavelets 

2.3.1  Introduction 

This  section  of  the  report  introduces  a  novel  approach  for  accomplishing  interactive  feature 
analysis  by  overcomplete  multiresolution  representations.  Traditional  wavelets  adapted  to 
“life  on  an  interval”,  can  overcome  “edge  effects”  of  wavelet  representations  on  a  line. 
Methods  of  contrast  enhancement  are  described  based  on  two  overcomplete  multiresolution 
representations  (interval  wavelets):  (1)  Deslanriers-Dubuc  interpolation,  (2)  Average 
interpolation. 
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We  show  quantitatively  that  transform  coefficients,  modified  by  an  adaptive  non-linear 
operator,  can  make  more  obvious  unseen  or  barely  seen  features  of  mammography  without 
requiring  additional  radiation. 

Arbitrary  regions  of  interest  (ROf)  of  a  mammogram  are  enhanced  by  average 
interpolation  (AI)  and  Deslauriers-Dubuc  interpolation  (DD)  representations  on  an 
interval.  The  results  of  local  (ROI)  enhancement  and  global  enhancement  of  mammograms 
are  compared  quantitatively.  We  demonstrate  that  our  method  can  provide  radiologists 
with  an  interactive  capability  to  support  high-speed  localized  processing  of  selected 
(suspicious)  areas  (lesions). 

We  demonstrate  that  features  extracted  from  multiscale  representations  can  provide  an 
adaptive  mechanism  for  accomplishing  local  contrast  enhancement.  By  improving  the 
visualization  of  breast  pathology  we  can  improve  chances  of  early  detection  while  requiring 
less  time  to  evaluate  mammograms  for  most  patients.  In  addition,  we  show  that  processing 
can  be  carried  out  at  real-time  (video  frame)  rates  over  selected  areas  of  arbitrary  shape. 
This  is  significant  in  consideration  of  tfie  large  image  matrix  size  of  digital  mammograms. 

2.3.2  Feature  analysis  via  interval  wavelets 

In  this  section,  we  describe  two  multiresolution  representations  which  were  investigated: 
average  interpolating  wavelet  transform  [31],  Deslauriers-Dubuc  interpolation  [32,  33]. 
These  representations  are  attractive  because  they  overcome  the  edge  effect  of  traditional 
multiresolution  representations  (based  on  periodization  of  a  finite  signal  to  a  signal  on  a 
line,  or  simply  adding  zeros  to  extend  a  signal  on  a  line,  etc).  The  shape  of  the  basis 
functions  for  these  representations  are  symmetric  or  antisymmetric,  and  allow  for  perfect 
reconstruction.  We  have  used  these  representations  to  decompose  an  arbitrary  rectangle  of 
a  mammogram,  so  that  the  rectangle  may  be  analyzed  independently. 

Cohen  and  Daubechies  [34]  first  adapted  multiresolution  analysis  on  the  line  to  “life  on 
the  interval” ,  where  a  sequence  of  successive  approximation  spaces  on  the  interval  were  be 
constructed  as:  \Jj  ezV  =  V([0,11).  n  ^.^2  Vj  =  {0}.  By  defining  Wj  as  an  orthogonal 
complement  of  Vj  in  V)_i,  Vj-i  =  Vj  0  Wj,  the  space  L^([0, 1])  can  be  represented  as  a 
direct  sum  L^([0, 1])  =  0jez  Wj. 

In  the  following  two  representations,  the  bases  are  compactly  supported,  but  not 
orthogonal.  These  representations  are  often  called  overcomplete  redundant  representations. 

Deslauriers-Dubuc  interpolation 

This  multiresolution  representation  consists  of  the  Deslauriers-Dubuc  fundamental 
functions  [32,  33].  Suppose  D  is  an  odd  integer  and  D  >  0,  a  fundamental  function  can  be 
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(a)  (b) 
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Figure  29:  (a)  Refinement  relation  for  Deslauriers-Dubuc  interpolation,  (b)  DD  wavelet  plot, 
D  =  3.  (c)  Refinement  relation  for  average  interpolation,  (d)  AI  wavelet  plot,D  =  2. 

defined  by  interpolating  the  Kronecker  sequence  recursively:  if  Fd  has  been  defined  in  the 
set  of  Bo  at  dyadic  rationals  fc/2°,  k  e  Z,  then  extend  Fd  to  the  set  of 

. .  ,Bj,j  e  Z.  Specifically,  when  the  value  of  the  function  is  defined  at  all  the 
dyadic  rationals  A:/2R  ^  G  Z,  we  extend  the  function  by  polynomial  interpolation  to  all 
dyadic  rationals  at  G  Z,  this  means  extend  the  function  to  all  the  dyadic  rationals 

half  way  between  the  previously  defined  positions.  When  j  oo^Fd  converges  to  a  unique 
continuous  function  on  the  real  line.  This  function  defines  an  (R,  D)  interpolating  wavelet 
for  R  =  R(T)),  and  is  compactly  supported.  Figure  29  shows  the  fundamental  solution  of 
DD  interpolation  and  associated  wavelets  {D  =  3). 

Donoho  [31]  showed  that  given  a  (R,  D)  interpolating  wavelet  tp,  we  can  construct  an 
interpolating  wavelet  transform,  mapping  the  function  /  into  approximation  sequence 
and  detail  sequence  djg+i^k,djg+2,k,  ■  ■  ■  ,dj^+rn.,k-,rn  G  Z,?-n  — )•  oo.  The  function  /  can 
then  be  reconstructed  from  its  coefficients, 

/  =  '^^^jo,k4’jo,k  +  ^  ^.i,k'>P.i,k- 

k  j>jo  k 

By  adapting  the  inhomogenous  interpolating  transform  to  “life  on  the  interval”,  we  can 
develop  additional  interpolating  wavelet  transforms  for  C[0, 1].  Suppose  that  is  a 
scaling  function  on  the  line.  The  scaling  functions  on  the  interval  are  derived  in  the 

following  way:  (1)  in  the  interior  of  the  interval,  they  are  just  the  same  as  on  the  line 

^  D  <k<2^  -D-l 
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(2)  on  the  edges  of  the  interval,  they  are  dilations  of  boundary  adjusted  functions: 


^nterv  ^  2^'/ -k),  0<k<  D, 


-  2^'  -  ^  -  1),  0  <  A;  <  D. 

Thus  for  the  spaces  Vj  [0, 1]  we  have  the  functions; 


inierv 


'  <5^Jf  {)<k<D 

'  4^j,k  D<k<2^  —  D  —  1 

fright  2i  _  D  _  1  <  A:  <  2T 


in  the  same  way,  we  can  define  the  wavelets  on  the  interval  for  the  detail  spaces  lTj[0, 1]: 


I  interv 
^hk 


'  <k  <  [D 1 2\ 

<  i>lk  \DI2\  <k  <2^  ~\DI2\ 
fright  2.  -lD/2\<k<  2T 


Thus  the  function  /  e  l'^j+i[0, 1]  can  be  written  as 


2^-1 


/  =  E  vi-e’"”  + 


k=0 


k=0 


interv 

k 


Average  interpolating  wavelet  transform 

Donoho  [31]  showed  that  connected  with  Deslauriers-Dubuc  interpolation,  there  exists 
an  average  interpolating  wavelet  transform.  Let  D  be  an  even  integer  and  >  0,  starting 
from  the  Kronecker  sequence  ao,k  =  ^k,0i  €  Z,  we  synthesize  the  averages  {aj,k)i.,j  £  Z+, 

the  sequence  Aj^d  =  E/t  cij,kl[k/2j,{k+-i)/2J]{t)  shall  converge  uniformly  to  a  fundamental 
solution  of  the  interpolating  scheme  Ad,  and  has  compact  support.  Figure  30  shows  the 
refinement  relation  and  its  associated  wavelets  with  D  =  2,  A. 

Aft)  may  be  written  as 

Aft)  =  J2^o,kADft  -  k). 

k 

Thus  the  average  interpolating  wavelet  transform  is  closely  related  with  Deslauriers-Dubuc 
interpolation.  Let  D  >  2  be  an  even  integer,  and  let  f  =  Ad  be  the  fundamental  solution 
of  the  average  interpolation,  and  =  Fd+\  (where  “DD”  represents  Deslauriers-Dubuc 
interpolation  )  be  a  fundamental  solution  of  the  Deslauriers-Dubuc  interpolation.  Let  tf  be 
the  AI  (  where  “AI”  represents  the  average  interpolation  )  wavelet,  and  let  be  the  DD 
wavelet.  Then  there  exists  a  relation; 

=  (i){x  +  1)  -  fix) 
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Mx)  i){x) 


D=4 

Figure  30:  Scaling  functions  and  associated  AI  wavelets. 
i,(x)  =  -  1/2)) 

as  shown  in  Figure  29.  If  the  (f)  and  ?/>  are  associated  with  the  average  interpolation,  we  can 
expand  the  function  /  €  Vji  into  a  sequence  of  coefficients  at  coarse  scale  in  the  spaces 
yjoijo  <  hi  9'i^cl  a  series  of  details: 

/  ”  ^jo,k4jo,k  d"  ^  ^ 

k  io<i<ii  k 

By  adapting  the  inhomogenous  interpolating  transform  on  an  interval”,  interpolating 
wavelet  transforms  for  L^[0, 1]  can  be  constructed.  For  the  spaces  V^  [0, 1]  we  have  the 
functions: 

'  0  <  fc  < 

^interv  ^  D<k<2^-D-1 

[fright  2J-F)-1<&<2F 

in  the  same  way,  we  can  define  the  wavelets  on  the  interval  for  the  detail  spaces  VFj[0, 1]: 

'  V’J/*  0<k<  D/2 

^xnterv  ^  D /2  <  k  <  2^  -  0/2 

[  i’it'  2^-Dl2<k<V. 

Thus  the  function  /  G  Vj+i[0, 1]  can  be  written  as 

/  =  E  +  E 

k=0  fc=0 

Figure  31  shows  plots  of  an  AI  refinement  relation  and  its  associated  wavelets  on  the 
interval. 
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(b) 


Figure  32:  (a)  Decomposition  and  enhancement  tree  structure  for  one  dimensional  case.  H 
and  G  are  a  filter  pairs.  E  is  a  nonlinear  operator  for  enhancement,  (b)  Selected  ROI  within 
a  mammogram,  (c)  ROI  is  processed  based  on  tensor  product:  each  row  is  processed  first, 
followed  by  the  processing  of  each  column. 
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Implementation 

Figure  32(a)  shows  the  one  dimensional  case  of  a  tree  structure  analysis  hlter  bank, 
used  to  implement  the  above  transformations.  After  decomposition  of  an  original  signal, 
enhancement  (dehned  by  E)  is  applied  to  wavelet  coefficients.  The  enhanced  signal  is  then 
obtained  by  reconstructing  the  signal  from  the  modified  coefficients.  By  extending  the 
processing  for  the  one  dimensional  case  to  two  dimensions,  we  were  able  to  enhance  features 
of  mammography  by  using  the  average  interpolation  and  Deslauriers-Dubuc  interpolation 
wavelet  transforms.  For  processing  of  an  arbitrary  region  in  a  mammogram,  we  used  a 
scanline  based  method:  First,  the  rows  of  a  selected  ROI  were  scanned  and  processed,  then 
the  colunms.  Figures  32(b)  and  (c)  illustrates  the  processing  steps  graphically. 

Enhancement  techniques 

To  accomplish  multiscale  contrast  enhancement,  non-linear  techniques  for  image 
enhancement  are  applied  to  the  multiresolution  representations.  In  our  case,  there  are  four 
components  in  the  transform  space;  horizontal,  vertical,  diagonal,  and  DC  component, 
represented  by  h\v\d\d  respectively,  where  i  is  the  level  of  a  transform.  Let  x  be  the 
original  mammogram,  /  be  the  function  designed  to  emphasize  features  of  importance 
within  a  selected  level  f,  T  be  the  total  level  of  transform.  Then  the  enhanced  image  may 
be  given  by 

i  =  (17) 

i=l 

In  general,  by  defining  function  g,  we  can  denote  specific  enhancement  schemes  for 
modifying  the  coefficients  within  distinct  levels  of  scale-space. 

There  are  three  ways  of  enhancement  techniques:  local  enhancement  technique  based 
on  multiscale  edges,  global  enhancement  techniques  of  multiscale  histogram  equalization 
and  multiscale  adaptive  gain  processing 

Experimental  results  and  discussion 

Preliminary  results  have  shown  that  the  multiscale  processing  techniques  described 
above  can  make  more  obvious  unseen  or  barely  seen  features  of  a  mammogram  without 
requiring  additional  radiation.  Our  study  suggests  that  the  analyzing  functions  presented 
in  this  paper  can  improve  the  visualization  of  features  of  importance  to  mammography  and 
assist  the  radiologist  in  the  early  detection  of  breast  cancer. 

Figure  33(a)  shows  a  “dense”  mammogram.  This  class  of  mammogram  is  more  typical 
in  younger  females  due  to  the  greater  absorption  of  X-ray  energy  by  less  fatty  tissues  in  the 
breast.  They  remain  particularly  difficult  to  diagnose  due  to  lack  of  contrast,  even  for 


48 


radiologists  specializing  in  mammography.  Figure  33(c)  shows  the  result  of  global  wavelet 
processing  for  four  levels  of  analysis.  In  this  case,  the  values  of  transform  coefhcients 
within  each  level  of  a  decomposition  (excluding  the  DC  cap)  were  modified  by  histogram 
equalization  independently.  Since  the  coefhcients  are  space-frequency  representations, 
contrast  modihcations  on  the  transform  side  are  preserved  in  part  on  the  spatial  side. 
Similar  contrast  gains  were  observed  for  additional  dense  radiographs.  Figure  33(b) 
displays  the  result  of  standard  histogram  equalization.  Unfortunately,  the  dense  tissues  of 
the  breast  image  are  “washed  out”  in  Figure  33(b). 

Mathematical  models  of  phantoms  were  constructed  to  validate  our  enhancement 
techniques  against  false  positives  arising  from  possible  artifacts  introduced  by  the  analyzing 
functions  and  to  compare  our  methods  against  traditional  image  processing  techniques  of 
improving  contrast.  Our  models  included  features  of  regular  and  irregular  shapes  and  sizes 
of  interest  in  mammographic  imaging,  such  as  microcalcihcations,  cylindrical  and  spicular 
objects,  and  conventional  masses.  Techniques  for  blending  a  normal  mammogram  with  the 
images  of  mathematical  models  were  developed.  The  purpose  of  these  experiments  was  to 
test  the  performance  of  our  processing  techniques  on  inputs  known  a  priori  using 
mammograms  where  the  objects  of  interest  were  deliberately  obscured  by  normal  breast 
tissues.  The  imaging  justification  for  blending  is  readily  apparent;  a  cancer  is  visible  in  a 
mammogram  because  of  its  (slightly)  higher  X-ray  attenuation  which  causes  a  lower 
radiation  exposure  on  the  film  in  the  appropriate  region  of  a  projected  image. 

Figure  34(b)  shows  an  example  of  a  mammogram  whereby  the  mathematical  phantom 
shown  in  Figure  34(c)  has  been  blended  into  a  clinically-proven,  cancer-free  mammogram. 
The  blended  image  was  shown  in  Figure  36(a),  it  was  constructed  by  adding  the  amplitude 
of  the  mathematical  phantom  image  in  Figure  34(c)  to  the  cancer  free  mammogram  in 
Figure  34(b)  followed  by  local  smoothing  of  the  combined  image. 

Radiologists  at  Shands  Hospital  at  the  University  of  Florida  validated  that  processing 
the  blended  mammogram  with  our  local  enhancement  techniques  introduced  no  signihcant 
artifacts  and  preserved  the  shapes  of  the  known  mammographic  features  (calcifications, 
dominant  masses,  and  spicular  lesions)  contained  in  the  original  mathematical  phantom. 

A  quantitative  measure  of  contrast  improvement  can  be  defined  by  a  Contrast 
Improvement  Index  (CII), 


CII  = 


C^Processed 


^  Ol’iginal 

where  Cprocessed  and  Coriginai  are  the  contrasts  for  a  region  of  interest  in  the  processed  and 
original  images,  respectively. 

In  this  study  we  adopt  a  version  of  the  optical  definition  of  contrast  introduced  by 
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Figure  33:  (a)  Original  dense  mammogram,  M56.  (b)  Enhancement  by  traditional  histogram 
equalization,  (c)  Global  enhancement  by  adaptive  histogram  equalization  of  Deslauriers- 
Dubuc  interpolation  (DD). 


(b)  (c) 


Figure  34:  (a)  Contrast  enhancement  by  adaptive  gain  processing  of  DD  interpolation 

wavelets,  (b)  Original  dense  mammogram,  M56.  (c)  Mathematical  phantom. 


Table  1:  Contrast  values  for  enhancement  by  local  enhancement  by  multiscale  edges  obtained 
from  average  interpolation  (AI),  Deslauriers-Dubuc  interpolation  (DD). 


Contrast  values  for  local  enhancement  techniques 

Feature 

C^Original 

Cai 

C'dd 

Minute  microcalcification  cluster 

0.0498 

0.1697 

0.1833 

Microcalcification  cluster 

0.0324 

0.0733 

0.1040 

Spicular  lesions 

0.0272 

0.0728 

0.0836 

Circular  (arterial)  calcification 

0.0378 

0.1060 

0.1242 

Well-circumscribed  mass 

0.0032 

0.0046 

0.0051 
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Figure  35:  Blended  mammogram:  (a)  Global  eiiliaiicemeiit  by  DD  gain,  (b)  Local  enliancement  by 
DD  edges,  (c)  Enhancement  of  rectangular  regions  by  DD  gain,  (d)  Enhancement  of  rectangular 
regions  by  DD  edge  processing. 


Table  2:  CII  for  enhancement  by  local  enhancement  by  multiscale  edges  obtained  from 
average  interpolation  (AI),Deslauriers-Dubuc  interpolation  (DD). 


Contrast  Improvement  Index  (CII)  for  local  enhancement  techniques 

Feature 

CIIai 

CIIdd 

Minute  microcalcification  cluster 

3.4051 

3.6773 

Microcalcification  cluster 

2.2605 

3.2061 

Spicular  lesions 

2.6756 

3.0728 

Circular  (arterial)  calcification 

2.8012 

3.2828 

Well-circumscribed  mass 

1.4364 

1.6024 

Morrow  et  al  [35].  The  contrast  C  of  an  object  is  defined  by 


where  /  is  the  mean  gray-level  value  of  a  particular  object  in  the  image,  called  the 
foreground,  and  b  is  the  mean  gray-level  value  of  a  surrounding  region  called  the 
background.  This  definition  of  contrast  has  the  advantage  of  being  independent  of  the 
actual  range  of  gray  levels  in  the  image.  With  the  aid  of  the  mathematical  phantom  we 
computed  local  masks  to  separate  the  foreground  and  background  regions  of  each  feature 
included  in  the  blended  mammogram. 

Figure  36(a)  shows  original  M56  (a  blended  mammogram).  Figure  35(a)  was  obtained 
by  global  enhancement  of  adaptive  gain  processing  of  Deslauriers-Dubuc  interpolation. 
Figure  34(a)  shows  enlarged  areas  containing  each  feature  in  the  processed  mammogram 
for  adaptive  gain  of  DD  interpolation  of  contrast  enhancement.  For  comparison  of  contrast, 
features  within  Figures  36(a)  and  35(a)  were  rescaled  collectively. 

Figure  35(a),  (b)  showed  enhancement  of  the  entire  mammogram,  and  Figure  35(c),  (d) 
showed  enhancement  of  the  five  pre-known  features  within  the  five  small  rectangles. 
Because  the  arbitrary  rectangles  enhancement  can  adaptively  select  parameters  according 
to  different  features,  it  is  more  flexible  than  the  enhancement  scheme  of  treating  all  the 
features  with  the  same  enhancement  parameters,  thus  enable  us  to  enhance  the  features 
more  effectively  while  reducing  the  enhancement  of  noise.  Figure  36  shows  the 
enhancement  of  an  arbitrary  region  of  interest  (ROI)  using  adaptive  gain  processing  of  DD 
interpolation  on  an  interval. 

By  constraining  the  enhancement  to  only  the  interest  region,  computation  is  greatly 
reduced.  Table  3  shows  the  comparison  of  actual  computation  time  of  processing  the  entire 
mammogram  vs  only  the  selected  ROI. 
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Figure  36:  Blended  mammogram:  (a)  Original  iiianimogram  blended  with  mathematical  phantom, 
i  b !  ROI  enhancement  by  adaptir'e  gain  processing  of  DD  inteiptolation  war'elet  transform. 


Table  3:  Comparison  of  computation  time.  TEntire-mammogram  represents  the  time  to  process 
a  complete  image  matrix,  while  Troi  represents  the  time  to  process  only  a  selected  ROI.  The 
number  of  pixels  within  the  ROI  shown  in  Figure  12  was  76,267  (executing  on  Sun  Sparc 
station  10/30). 


Computation  time  (in  seconds)  comparison  of  whole  mammogram  vs  ROI 

Matrix  size  (number  of  pixels) 

TEntire-mammogram 

Troi 

TEntire-mammogram  1  TflO  I 

512x512 

748 

135 

5.54 

1024x1024 

5760 

135 

42.67 

2.3.3  Auto-correlation  shell  representations 

This  section  introduces  another  approach  for  accomplishing  mammographic  feature 
analysis  by  an  overcomplete  multiresolution  representation:  auto-correlation  shell 
representation  of  compactly  supported  wavelets.  We  show  that  an  arbitrary  region  of 
interest  (ROI)  of  the  mammogram  can  be  enhanced  using  the  auto-correlation  shell.  The 
ROI  enhancement  provides  the  radiologists  a  way  to  view  the  mammogram  with  only  the 
ROI  enhanced  on  the  screen. 

Analytical  formulation  for  auto-correlation  shells 

In  this  section,  we  have  investigated  auto-correlation  functions  of  Daubechies’s 
compactly  supported  wavelets  [36].  This  representation  is  attractive  because  the  shape  of 
the  basis  functions  for  this  representation  are  symmetric,  and  allow  for  perfect 
reconstruction.  We  have  used  this  representation  to  decompose  an  arbitrary  region  of 
interest  (ROI)  within  a  mammogram,  so  that  the  ROI  may  be  analyzed  independently. 

Since  the  coefficients  of  orthogonal  wavelet  expansions  are  not  shift  invariant  in  general, 
it  can  be  quite  difficult  to  explore  the  property  of  features  from  scale  to  scale.  The 
asymmetric  shape  of  compactly  supported  wavelets,  makes  these  wavelets  poor  edge 
detectors,  and  can  introduce  artifacts  when  image  enhancement  is  based  on  wavelet 
coefficients  alone.  However  the  compact  support  characteristic  of  these  bases  make  them 
exact  both  in  decomposition  and  reconstruction.  The  auto-correlation  shell  basis  function 
has  the  same  support  width  of  its  associated  compactly  supported  wavelets. 

To  overcome  the  drawbacks  of  compactly  supported  orthogonal  wavelets  in  the 
application  of  image  analysis,  Saito  and  Beylkin  [36]  introduced  the  notion  of  an 
auto-correlation  shell  .  Let  ip{x)  be  some  compactly  supported  wavelet,  and  (^{x)  be  the 
scaling  function.  By  dehnition  of  auto-correlation,  we  can  write 

/4-00  r+oo 

Hy)Hy  -  ^)dy,  T'(x)=/  xl;{y)xp{y  -  x)dy 

-oo  J  — oo 
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where  their  respective  Fourier  transforms  are 

^(0  =  \m\\  m  =  mw- 

Both  and  v];(a;)  are  supported  within  the  interval  [—L  +  1,  L  +  1],  where  L  is  the 
length  of  the  quadrature  mirror  filter  of  a  compactly  supported  orthogonal  wavelet,  and 
vl/(a:;)  and  ^>{x)  are  symmetric.  Figure  38  and  Figure  39  show  the  analyzing  filters  of  the 
auto-correlation  shell  and  Daubechies’s  wavelet  with  two  vanishing  moments. 

The  functions  ^(x)  may  be  viewed  as  a  pseudo-differential  operator,  this  allow  us  to 
locate  the  edges  of  a  signal  at  different  scales  of  auto-correlation  expansions.  From  Figure 
40,  we  can  see  that  the  auto- correlation  functions  #(x)  and  ^(x)  are  smoother  than  the 
functions  <?!)(x)  and  V^(x),  and  and  ’F(^)  decay  faster  than  and  respectively. 

The  auto-correlated  shell  basis  function  of  Daubechies’s  compactly  supported  wavelets  with 
N  =  2  are  shown  in  figure  37. 

Let  {hk}-L+i<k<L-i  be  the  low  pass  hlter  related  with  $(x),  and  let  {gk}-L+i<k<L-i  be 
the  high  pass  filter  related  to  T(x).  We  have  implemented  a  pyramid  algorithm  for 
expanding  a  function  /  into  an  auto-correlation  shell: 

L-l 

l=-L+l 

L-l 

d'j,k  =  9ifj-iM2}-'^i- 

l=-L+l 

Reconstruction  is  accomplished  by 

fj-l,k  — 


2.3.4  Summary 

In  this  study,  methods  for  accomplishing  adaptive  contrast  enhancement  by  multiscale 
representations  were  further  investigated.  Contrast  enhancement  was  applied  to  features  of 
specific  interest  to  mammography  including  masses,  spicules  and  microcalcifications. 
Multiresolution  representations  provided  an  adaptive  mechanism  for  the  local  emphasis  of 
such  features  blended  into  digitized  mammograms.  Average  interpolation  and 
Deslauriers-Dubuc  interpolation  representations  on  an  interval  enabled  us  to  enhance 
arbitrary  regions  of  interest  within  a  mammogram.  The  enhancement  of  ROFs  provides  the 
radiologists  an  interactive  way  of  processing  only  an  interesting  suspicious  region  of  a 
mammogram,  while  reducing  the  computational  cost  compared  to  processing  an  entire 
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Figure  37:  Plots  of  auto-correlation  functions  of  Daubechies’s  compactly  supported  wavelet 
with  b  =  4.  (a)  $(x)  ,  (b) 
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Figure  38:  Analyzing  filters  for  auto-correlation  shell  of  Daubechies’s  compactly  support 
wavelets  with  two  vanishing  moments  (three  levels  of  analysis  shown). 
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Figure  39:  Analyzing  filters  for  Daubechies’s  compactly  supported  wavelets  with  two  van¬ 
ishing  moments  (three  levels  of  analysis  shown). 


Figure  40:  1-D  signal  expanded  with  an  auto-correlation  shell.  The  locations  of  the 

edges  in  the  original  signal  So  correspond  to  the  zero-crossings  in  the  detail  information 
Di.  D2,  Ds-,  D4  . 


mammogram.  In  addition,  the  representation  of  an  auto- correlation  shell  was  studied.  This 
representation  can  also  be  used  to  enhance  mammographic  features  over  an  arbitrary 
region.  These  initial  results  are  encouraging  and  suggest  that  wavelet  based  image 
processing  algorithms  could  play  an  important  role  in  improving  the  imaging  performance 
of  digital  mammography. 

2.4  Quantitative  Evaluation  of  Clinical  Images 

2.4.1  Introduction 

This  part  of  the  report  relates  to  the  development  of  objective  ways  to  assess  the 
performance  of  wavelet  image  processing  algorithms.  The  objective  is  to  develop  techniques 
to  evaluate  wavelet  algorithms  so  they  can  then  be  optimized  for  clinical  use  in 
m  ammogr  ap  hy. 

Substantial  progress  has  been  made  in  the  development  of  techniques  to  assist  in  the 
quantitative  evaluation  of  wavelet  based  image  processing  algorithms.  In  addition,  these 
techniques  have  been  applied  to  optimize  the  three  parameters  available  in  current  wavelet 
based  algorithms.  Specific  achievements  in  the  past  year  include; 

•  a  demonstration  of  the  ability  of  optimized  wavelet  based  algorithms  to  make  visible 
simple  objects  in  a  noisy  background  which  were  previously  invisible; 

•  a  demonstration  of  the  inherent  superiority  of  wavelet  based  algorithms  for  the 
detection  of  simple  objects  as  compared  to  algorithms  frequently  used  in  medical 
imaging  including  unsharp  mask  enhancement  and  median  filtering; 

•  a  demonstration  of  the  special  requirements  of  wavelet  algorithms  when  enhancing 
the  visibility  of  features  of  specific  interest  in  mammography,  namely 
microcalcifications,  masses  and  fibril  structures. 

2.4.2  Specific  projects  completed 

1.  Image  quality  index.  Computer  simulated  phantoms  are  an  attractive  choice  for 

evaluation  of  image  processing  algorithms  since  the  features  of  interest  are  known  and 
well  quantihed.  We  developed  a  metric  of  the  improvement  in  image  Signal  to  Noise 
Ratio  (SNR)  and  used  an  Enhancement  Factor  (EF  =  SNRo/SNRi)  where  the 
subscripts  o  and  i  refer  to  the  processed  image  and  original  image,  respectively.  We 
validated  that  there  was  an  excellent  correlation  between  the  computed  EF  and 
observer  performance  in  psychophysics  tests  (see  reference  1  in  section  2.5). 
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2.  Optimization  of  wavelet  algorithms.  Three  separate  wavelet  basis  functions  were 
investigated  (i.e.  dyadic,  hexagonal  and  phi).  Each  algorithm  has  three  parameters 
(Level  [L],  Threshold  [T],  and  Gain  [G])  which  may  be  varied.  The  effects  of  these 
three  parameters  on  all  three  wavelet  basis  functions  were  systematically  investigated 
for  simple  structures  (e.g.  Gaussian  signals)  embedded  in  white  (Gaussian)  noise. 

This  work  concluded  that  the  three  parameters  have  similar  effects  irrespective  of  the 
selected  wavelet  basis.  In  addition,  the  study  demonstrated  that  wavelet  based 
algorithms  are  generally  not  linear  and  that  significant  improvements  in  SNR  can  be 
achieved  with  optimized  values  of  the  parameters  L,  T  and  G.  Another  encouraging 
Ending  was  that  the  false  positive  rate  was  generally  very  low  provided  that  the  gain 
parameter  was  not  set  too  high  (see  references  1,  3,  6  and  9  in  section  2.5). 

3.  Size  and  shape  factors.  Two  important  variables  for  the  objects  of  interest  in 
mammograms  are  object  shape  and  object  size.  This  is  of  particular  relevance  given 
the  inherent  multiscale  capabilities  of  wavelet  signal  decomposition. 

The  importance  of  shape  was  systematically  investigated  and  published  (see  reference 
3  in  section  2.5)  and  a  paper  describing  in  detail  the  significance  of  shape  will  be 
presented  at  the  forthcoming  SPIE  meeting  in  San  Diego  in  Eebruary  1995  (see 
reference  6  in  section  2.5).  In  general,  the  size  of  a  feature  of  interest  is  more 
important  than  its  specific  shape.  Eurthermore,  it  is  possible  to  optimize  wavelet 
algorithms  for  any  selected  feature  size. 

4.  Mammographic  features.  The  specific  objects  of  interest  in  mammograms  are  masses, 
microcalcification  clusters  and  fibrils.  The  ability  of  wavelet  based  algorithms  to 
enhance  these  features  were  investigated  using  both  computer  simulated  phantoms, 
and  physical  phantoms  radiographed  using  a  screen/film  combination  designed  for 
mammography  and  subsequently  digitized  to  permit  image  processing. 

Results  of  these  investigations  have  been  published  for  both  phantom  images  (see 
references  4  and  9)  and  for  radiographs  of  the  American  College  of  Radiology  (ACR) 
accreditation  phantom  (see  reference  2  in  section  2.5).  These  studies  showed  that  the 
differing  types  of  features  investigated  would  benefit  significantly  from  specially 
tailored  algorithms. 

5.  Clinical  applications.  In  addition  to  the  phantom  studies,  we  have  performed 
investigations  into  the  effectiveness  of  digitized  clinical  mammograms.  These  studies 
have  generally  shown  that  the  trends  identified  with  computer  simulated  physical 
phantoms  and  with  radiographs  of  the  ACR  phantom  are  similar  to  those  found  with 
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clinical  mammograms.  The  most  promising  wavelet  basis  identified  by  all  three 
approaches  is  the  dyadic  basis.  The  enhancement  strategy  that  shows  the  most 
promise  is  an  Edge  enhancement  where  it  is  only  the  wavelet  maxima  at  a  selected 
scale  which  are  enhanced  by  a  gain  factor.  Two  alternative  wavelet  coefficient 
enhancement  strategies  investigated  (Gain  and  Histogram  Equalization)  were 
generally  found  to  be  inferior  to  Edge  enhancement.  Full  details  of  these  results  are 
available  in  the  published  papers  (see  references  2,  7  and  8). 

2.4.3  Current  projects 

1.  Image  display.  We  are  currently  investigating  ways  of  presenting  digital 
mammographic  images  on  soft  copy  displays.  This  issue  is  difficult  because  of  the 
large  matrix  sizes  and  the  number  of  bits  used  to  code  for  any  individual  pixel.  Our 
preliminary  results  were  presented  at  the  SPIE  meeting  on  Medical  Imaging  in  San 
Diego  in  February  1994  (see  reference  5  in  section  2.5). 

2.  Digital  x-ray  images  (LoRad  system).  A  recently  installed  LoRad  biopsy  device 
incorporates  a  small  field  of  view  high  resolution  (10242)  digital  imaging  device.  The 
clinical  use  of  this  system  is  to  permit  needle  core  biopsies  to  be  obtained.  We  are 
currently  using  this  system  to  radiograph  images  of  the  ACR  phantom  and  also  using 
clinical  images,  (see  reference  5  in  section  2.5). 

3.  Structured  background.  There  are  two  major  sources  of  noise  in  mammograms.  The 
first  arises  from  quantum  mottle  and  film  granularity,  and  this  type  of  noise  may  be 
simulated  using  the  computer  generated  phantom  images  as  described  above.  A 
second  source  of  ’’noise”  is  the  structure  patterns  in  mammograms.  A  major  thrust  of 
the  work  being  planned  is  to  investigate  how  a  structured  mammographic  background 
affects  the  visibility  of  features  which  identify  cancers  (masses,  microcalcifications 
and  fibrils).  In  addition,  the  ability  of  wavelet  based  algorithms  to  improve  feature 
visibility  will  be  investigated.  Of  particular  importance  is  whether  algorithms 
optimized  for  random  Gaussian  noise  are  also  optimized  for  structured  background. 

4.  Digital  mammogram  databases 

We  have  purchased  a  Mammographic  database  from  the  Mammographic  Image 
Analysis  Society  in  England  with  consists  of  322  images  digitized  with  a  50 
micrometer  spots  size.  In  addition,  we  are  in  the  process  of  digitizing  a  collection  of 
biopsy  proven  cases  from  the  University  of  Florida  to  use  for  evaluation  in  our  studies 
of  wavelet  analysis. 
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5.  Wavelet  compression 

We  are  in  the  process  of  developing  methods  for  applying  wavelet  and  jpeg 
compression  to  mammograms  with  the  intention  of  evaluating  both  methods  for  the 
best  compression  method  for  these  large  image  files.  Receiver  Operator 
Characteristic  (ROC)  analysis  as  well  as  forced  choice  analysis  will  be  carried  out  to 
determine  the  best  method  and  the  limits  of  compression. 

6.  Sampling  images  for  display 

We  are  beginning  the  investigation  of  subsampling  methods  for  reducing  image 
matrix  size  for  digital  mammographic  images.  We  will  evaluate  methods  using 
Receiver  Operating  Characteristics  and  forced  choice  analysis  as  well  as  physical 
measurements  to  determine  the  best  clinical  way  to  subsample  as  well  as  to  develop 
methods  for  correlating  physical  measurements  with  clinical  accuracy. 

7.  Mammography  database 

We  are  updating  the  mammography  database  to  meet  the  American  College  of 
Radiology’s  new  rules  for  outcomes  reporting. 

2.5  Published  journal  papers  and  conference  proceedings 

1.  Xing  Y,  Huda  W,  Laine  A,  Fan  J  ’’Simulated  phantom  images  for  optimizing  wavelet 
based  image  processing  algorithms  in  mammography”  Proceedings  of  the  SPIE,  San 
Diego,  25-26  July  1994,  Volume  2299  207-217. 

2.  Qu  G,  Huda  W,  Laine  A,  Steinbach  BG,  Honeyman  JC  ’’Use  of  accreditation 
phantoms  and  clinical  images  to  evaluate  mammography  image  processing 
algorithms”  In  Digital  Mammography,  Editors  AG  Gale  et  al,  Elsevier  Science  BV; 
Amsterdam  (1994)  345-  354. 

3.  Jing  Z,  Zheng  Y,  Huda  W,  Laine  A  and  Fan  J  ’’Mathematical  models  for  quantitative 
evaluation  of  wavelet  based  image  processing  algorithms”  Proceedings  of  the  SPIE, 
San  Diego,  25-26  July  1994,  SPIE  Volume  2303  569-578. 

4.  Laine  A,  Schuler  S,  Fan  J  and  Huda  W  ’’Mammographic  feature  enhancement  by 
multiscale  analysis”,  IEEE  Transactions  on  Medical  Imaging,  Vol.  13,  No. 4, 
December,  1994. 
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5.  Honeynian  JC,  Chen  D,  Huda  W  and  Steinbach  BG  ’’Impact  of  image  display 
parameters  on  the  visibility  of  microcalcifications  and  masses”  Submitted  to  the 
SPIE. 

6.  Huda  W,  Xing  Y,  Lain  A,  Fan  J,  Steinbach  BG  and  Honeyman  JC  ’’Assessment  of  a 
wavelet  image  processing  algorithm  for  mammography”  Submitted  to  SPIE. 

7.  Steinbach  BG,  Laine  A,  Huda  W  and  Honeyman  JC  ’’Mammography  image 
enhancement  by  using  wavelet  transforms”  Submitted  to  RadioGraphics. 

8.  Laine  A,  Huda  W,  Steinbach  BG  and  Honeyman  JC  ’’Mammographic  image 
processing  using  wavelet  processing  techniques”  submitted  to  European  Radiology. 

9.  Xing  Y,  Huda  W,  Laine  A  and  Fan  J  ’’Optimization  of  a  dyadic  wavelet  for 
mammography”  In  preparation  for  submission  to  IEEE  Transactions  on  Medical 
Imaging. 

2.6  Published  abstracts  and  presentations 

1.  Title:  Use  of  phantoms  and  clinical  images  to  evaluate  wavelet  based  image 
processing  algorithms  in  mammography. 

Authors:  G  Qu,  W  Huda,  A  Laine,  BG  Steinbach  and  JC  Honeyman. 

Presented  at  the  1994  American  Association  of  Physicists  in  Medicine  (AAPM) 
meeting  in  Anaheim,  CA 

Abstract:  Wavelet  based  image  processing  algorithms  have  been  proposed  to  enhance 
clinically  important  features  in  mammograms.  In  this  study,  images  of  an  ACR 
mammographic  phantom  and  clinical  images  were  both  used  to  evaluate  the 
subjective  improvement  of  image  quality.  Readers  were  asked  to  compare  the 
visibility  of  masses,  spicules  and  microcalcifications  of  images  processed  using  nine 
different  wavelet  algorithms.  The  visibility  rankings  ranged  from  a  value  of  1 
(markedly  less  than  the  original)  to  5  (markedly  better  than  the  original).  For  the 
ACR  phantom  images,  a  wide  range  of  positive  and  negative  rankings  were  obtained 
suggesting  that  algorithms  could  be  very  good  or  very  bad.  For  clinical  images, 
however,  most  rankings  tended  to  correspond  to  the  mean  value.  The  advantages  and 
limitations  of  each  phantom  will  be  discussed. 
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2.  Title:  Investigation  of  wavelet  based  image  processing  algorithms  for  use  in 
mammography 

Authors:  Y  Zheng,  Y  Xing,  W  Huda  and  A  Laine. 

Presented  at  the  1994  American  Association  of  Physicists  in  Medicine  (AAPM) 
meeting  in  Anaheim,  CA 

Abstract:  It  is  important  to  optimize  image  processing  algorithms  prior  to  subjecting 
them  to  clinical  evaluation.  In  this  study,  computer  simulated  phantoms  were  used  to 
systematically  investigate  the  parameters  of  wavelet  image  processing  algorithms. 
These  included  the  hierarchical  decomposition  level,  wavelet  coefficient  thresholds 
and  enhancement  gain  values  [Laine  et  al,  SPIE  Volume  1905  (1993)  521-532)]. 

Image  processing  algorithms  designed  for  use  in  mammography  were  applied  to  these 
computer  generated  images  and  the  input/output  SNR  values  determined  as  a 
function  of  each  parameter.  This  methodology  permits  the  identification  of  promising 
wavelet  parameters  which  could  permit  mammograms  to  be  processed  to  improve  the 
visibility  of  features  such  as  masses,  spiculations  and  microcalcifications. 

3.  Title:  Assessment  of  contrast  enhancement  for  wavelet  based  mammography  image 
processing  algorithms. 

Authors:  W  Huda,  A  Laine,  JC  Honeyman,  BG  Steinbach. 

Presented  at  the  1994  Canadian  Organization  of  Medical  Physics  (COMP)  meeting  in 
Toronto 

Abstract:  Wavelet  based  image  enhancement  algorithms  have  been  proposed  for  use 
in  mammography  to  improve  visibility  of  clinically  important  features.  In  this  study, 
three  wavelet  based  algorithms  were  evaluated  for  their  ability  to  enhance  the 
visibility  of  masses,  spicular  lesions  and  microcalcifications  which  had  been 
mathematically  blended  into  digitized  clinical  mammograms.  The  results  were 
compared  with  a  competitive  local  enhancement  (unsharp  masking)  algorithm. 
Contrast  Improvement  Indices  (Oils)  were  obtained  by  measuring  the  ratio  of 
contrast  in  the  processed  image  to  that  of  the  original  image.  For  the  wavelet 
algorithms  investigated,  CII  values  were  in  the  range  1.9  to  4.5  for  the  three  types  of 
added  mammographic  feature.  The  corresponding  CII  values  for  the  unsharp  mask 
algorithm,  however,  were  in  the  range  1.1  to  1.8.  The  results  of  this  study  indicate 
that  wavelet  based  image  processing  algorithms  could  be  useful  for  enhancing  clinical 
mammograms . 
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4.  Title:  Mammography  image  enhancement  using  wavelet  transforms. 

Authors:  BG  Steinbach,  A  Laine,  JC  Honeyman,  W  Huda. 

Presented  at  the  meeting  of  the  American  Rontgen  Ray  Society 

Abstract:  We  present  a  novel  approach  for  accomplishing  mammographic  feature 
analysis  utilizing  wavelet  transforms.  Identification  of  image  features  at  different 
scales  permits  selected  attributes  in  a  mammogram  to  be  enhanced.  Using  digitized 
mammograms,  wavelet  transform  algorithms  were  evaluated  in  several  cases  of  biopsy 
proven  breast  carcinomas.  Spiculated  and  well  defined  masses  of  varying  degrees  of 
subtlety  were  processed  using  image  processing  algorithms  based  on  multiscale 
wavelet  transforms.  For  comparison  purposes,  these  images  were  also  processed  using 
traditional  histogram  equalization. 

Mammograms  processed  with  algorithms  based  on  the  wavelet  transform  showed 
some  improvements  in  feature  visualization.  This  Wavelet  Transform  method  has  the 
potential  to  significantly  improve  the  visibility  of  clinically  important  features  (i.e. 
masses,  spicules  &  microcalcifications)  and  merits  further  investigation. 

5.  Title:  Assessment  of  a  wavelet  image  processing  algorithm  for  mammography. 

Authors:  W  Huda,  Y  Xing,  A  Laine,  J  Fan,  BG  Steinbach,  JC  Honeyman. 

Submitted  to:  Image  Processing,  Murray  H  Loew,  at  Ml 95  Medical  Imaging  1995,  26 
February  -  2  March  1995,  San  Diego 

Purpose:  Image  processing  algorithms  based  on  the  wavelet  transform  has  been 
proposed  to  improve  the  visibility  of  mammographic  features.  In  this  study,  a 
computer  simulated  phantom  was  used  to  optimize  a  dyadic  wavelet  algorithm.  The 
degree  of  image  enhancement  achieved  with  this  algorithm  was  compared  with 
algorithms  currently  used  in  diagnostic  radiology. 

Method:  A  computer  simulated  phantom  was  generated  which  contained  masses, 
fibrils  and  microcalcifications  together  with  added  random  noise.  Image  improvement 
was  determined  from  the  ratio  of  output  to  input  signal  to  noise  ratios.  Performance 
of  the  following  three  algorithms  was  investigated:  1)  dyadic  wavelet  transform;  2) 
unsharp  masking;  3)  median  filtering. 

Results:  Each  algorithm  has  several  free  parameters  which  may  be  altered  and  affect 
the  degree  of  image  enhancement  achieved.  The  properties  of  these  parameters  for 
each  algorithm  were  studied.  Optimal  parameters  for  these  algorithms  were  found  to 
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depend  on  the  type  of  mammographic  feature  being  studied. 

Conclusions:  The  results  obtained  show  that  a  multiresolution  approach  to 
mammographic  image  processing  using  a  dyadic  wavelet  transform  is  generally  a 
more  powerful  tool  than  those  offered  by  traditional  image  processing  algorithms. 

6.  Title:  Impact  of  image  display  parameters  on  the  visibility  of  microcalcifications  and 
masses. 

Authors:  JC  Honeyman,  D  Chen,  W  Huda,  BG  Steinbach. 

Submitted  to:  Image  Dis-plaij,  Yongmin  Kim,  at  Ml 95  Medical  Imaging  1995,  26 
February  -  2  March  1995,  San  Diego 

Purpose:  High  resolution  digital  mammography  systems  are  likely  to  appear  in 
clinical  practice  in  the  near  future.  Digital  image  data  will  permit  the  manner  in 
which  mammograms  are  displayed  to  be  readily  modihed.  In  this  study,  the 
importance  of  display  parameters  on  mammography  feature  conspicuity  was 
evaluated  using  images  printed  on  film  and  displayed  on  diagnostic  workstations. 

Method:  Clinical  mammographic  images,  5  x  5  cm  in  size,  were  obtained  from  a 
LoRad  Digital  Spot  Mammography  system.  The  image  matrix  size  was  5122  with  a 
pixel  dimension  of  about  100  em.  Ten  images  had  microcalcifications  and  ten  images 
had  masses.  For  each  image,  nine  copies  were  made  with  differing  window  and  level 
settings  which  were  printed  on  film  and  displayed  on  a  monitor.  Five  radiologists 
specializing  in  mammography  ranked  conspicuity  of  each  mammographic  feature  from 
the  least  visible  to  the  most  visible  in  a  9  forced  choice  experiment. 

Results:  The  characteristic  curves  for  transforming  image  pixel  intensity  to  optical 
density  (film)  and  image  monitor  intensity  (workstations)  were  obtained.  The 
significance  of  window  and  level  on  the  perceived  visibility  of  mammographic  features 
were  determined.  There  are  differences  between  small  high  contrast  features  such  as 
mammographic  calcifications  and  larger  low  contrast  features  such  as  malignant 
masses.  The  relative  ranking  of  the  conspicuity  of  mammographic  features,  however, 
were  generally  similar  for  both  film  and  workstation  displays. 

Conclusions:  The  data  obtained  quantify  the  significance  of  image  display  parameters 
for  improving  the  conspicuity  of  mammographic  features.  These  results  will  permit 
the  optimization  for  printing  films  from  digital  mammography  data  and  assist  in  the 
design  of  display  workstations  for  use  in  mammography. 
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2.7  Lectures  and  invited  talks 

1 .  Andrew  Laine,  International  Conference  of  the  IEEE  Engineering  in  Medicine  and 
Biology  Society  (EMBS),  November  1-2,  1994,  Baltimore,  MD,  Invited  speaker. 

2.  Andrew  Laine,  Yale  University,  Department  of  Mathematics,  Colloquium  guest 
speaker,  April,  1994. 

3.  Andrew  Laine,  University  of  Chicago,  Departments  of  Computer  Science  and 
Radiology,  Invited  lectures.  Colloquium  speaker,  January,  1994. 
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3  Conclusions 


During  the  past  year,  we  have  made  significant  progress  in  the  development  of  a 
methodology  for  accomplishing  adaptive  contrast  enhancement  by  multiscale 
representations.  Our  studies  have  demonstrated  that  features  extracted  from 
multiresolution  representations  can  provide  an  adaptive  mechanism  for  the  local  emphasis 
of  salient  and  subtle  features  of  importance  to  mammography.  The  improved  contrast  of 
mammographic  features  make  these  techniques  apijealing  for  computed  aided  diagnosis 
(CAD)  and  screening  mammography.  Screening  mammography  examinations  are  certain  to 
grow  substantially  in  the  next  few  years,  and  analytic  methods  that  can  assist  general 
radiologists  in  reading  mammograms  shall  be  of  great  importance. 

In  the  paragraphs  below,  we  summarize  our  progress  and  identify  future  directions  of 
research  to  be  carried  out  during  the  next  year  of  our  investigation. 

We  have  established  connections  between  dyadic  wavelet  enhancement  algorithms  and 
traditional  unsharp  masking.  We  proved  that  two  cases  of  linear  enhancement  were 
mathematically  equivalent  to  traditional  unsharp  masking  with  Gaussian  low-pass  filtering. 
We  designed  a  methodology  for  nonlinear  enhancement  with  a  simple  nonlinear  function  to 
overcome  the  dynamic  range  requirement  usually  associated  contrast  enhancement  of 
digital  radiographs.  By  careful  selection  of  wavelet  filters  and  enhancement  function,  we 
showed  that  artifacts  can  be  eliminated.  An  additional  advantage  of  our  simple 
enhancement  function  is  that  it  includes  traditional  unsharp  masking  as  a  subset.  We 
showed  how  an  edge-preserved  denoising  stage  (wavelet  shrinkage)  can  be  appropriately 
incorporated  into  our  contrast  enhancement  framework,  and  introduced  a  method  for 
adaptive  threshold  estimation.  We  showed  how  denoising  and  enhancement  operations 
should  be  carried  out  for  two  dimensional  images  to  avoid  orientation  distortions. 

Our  future  research  plans  include  the  systematic  study  of  gain  and  threshold  parameters 
for  the  nonlinear  enhancement.  In  addition,  in  the  next  year  we  shall  seek  localized  and 
complex  nonlinear  methods  to  improve  the  performance  of  our  existing  algorithm. 

We  have  design  and  implemented  an  interactive  image  editor  for  digital  mammography, 
called  ”Xmam”.  We  have  defined  the  functionality  and  data  structures  of  a  computer 
image  editor  that  can  allow  radiologists  to  interactively  indicate  on  a  computer  screens 
(simultaneously  displaying  four  radiographic  views  of  screening)  regions  diagnosed  as 
probable  cancer.  This  tool  shall  facilitate  the  compilation  of  a  large  database  of  diagnosed 
cases  to  provide  ground  truth  for  the  development  and  verification  of  our  computer 
algorithms  for  the  detection  of  breast  cancer. 
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During  the  next  year,  XMam  shall  evolve  into  a  front-end  for  our  state-of-the-art 
mammography  workstation  as  additional  software  is  developed  and  we  migrate  our  system 
onto  more  sophisticated  hardware. 

In  addition,  methods  for  accomplishing  adaptive  contrast  enhancement  by  multiscale 
representations  were  further  investigated.  Contrast  enhancement  was  applied  to  features  of 
importance  to  mammography  including  masses,  spicnles  and  microcalcifications. 
Multiresolution  representations  provided  an  adaptive  mechanism  for  the  local  emphasis  of 
such  features  blended  into  digitized  mammograms.  Average  interpolation  and 
Deslauriers-Dubuc  interpolation  representations  on  an  interval  enabled  us  to  enhance 
arbitrary  regions  of  interest  within  digital  mammograms. 

Enhancement  of  arbitrarily  shaped  ROPs  provides  the  radiologist  with  an  interactive 
capability  to  process  only  suspicious  regions  of  a  mammogram,  while  significantly  reducing 
execution  time  compared  to  processing  an  entire  image  matrix.  In  addition,  the 
representation  of  an  auto-correlation  shell  was  stndied.  We  remain  excited  about  this 
representation  and  plan  to  further  study  its  ability  to  enhance  mammographic  features 
over  arbitrary  regions.  These  initial  results  are  encouraging  and  suggest  that  wavelet  based 
image  processing  algorithms  shall  play  an  important  role  in  improving  the  imaging 
performance  of  digital  mammography. 

Finally,  we  reported  on  the  development  of  objective  ways  to  assess  the  performance  of 
wavelet  image  processing  algorithms.  Our  objective  is  to  develop  techniques  to  evaluate 
wavelet  algorithms  so  they  can  then  be  optimized  for  clinical  use  in  mammography. 

Substantial  progress  was  made  in  the  development  of  techniques  to  assist  in  the 
quantitative  evaluation  of  wavelet  based  image  processing  algorithms.  In  addition,  these 
techniques  were  applied  to  optimize  the  three  parameters  available  in  current  wavelet  based 
algorithms. 

During  the  next  year,  we  expect  to  continue  our  study  of  wavelet  filter  design  and 
selection  in  the  rehnement  of  our  early  basis  functions  (dyadic  wavelets,  phi-transform,  and 
hexagonal  wavelets)  conceived  as  initial  instances  in  the  evolution  of  three  specialized 
detectors.  As  described  in  our  research  plan,  we  shall  design  three  ’’detectors”  to  focus  on 
three  distinct  types  of  mammographic  features:  (1)  microcalcihcations,  (2)  spicular  lesions 
and  (3)  masses. 

In  summary,  we  have  exceeded  the  goals  as  described  in  our  Statement  of  Work  for  the 
second  year,  and  our  research  and  development  plans  remain  on  schedule. 
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